######################################
## 18_reportPartyPressure_JOP_RnR.R ##
######################################

rm (list=ls())

print (version)

# platform       aarch64-apple-darwin20      
# arch           aarch64                     
# os             darwin20                    
# system         aarch64, darwin20           
# status                                     
# major          4                           
# minor          1.2                         
# year           2021                        
# month          11                          
# day            01                          
# svn rev        81115                       
# language       R                           
# version.string R version 4.1.2 (2021-11-01)
# nickname       Bird Hippie 

library (xtable)
library (sp)
library (coda)


###############################
#### Brazil party pressure ####
###############################
rm (list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd() 

# Load all party pressure objects
load (file="brazil_ptirt_mcmc.Rda")

# Check convergence
gD <- gelman.diag (chainsPTIRT, multivariate=F) 
summary (gD[[1]][,1])

Beta.e   <- rbind (chainsPTIRT[[1]][,grep("beta",  colnames (chainsPTIRT[[1]]))],
                   chainsPTIRT[[2]][,grep("beta",  colnames (chainsPTIRT[[2]]))])

Alpha.e   <- rbind (chainsPTIRT[[1]][,grep("alpha",  colnames (chainsPTIRT[[1]]))],
                    chainsPTIRT[[2]][,grep("alpha",  colnames (chainsPTIRT[[2]]))])

X.e   <- rbind (chainsPTIRT[[1]][,grep("X",  colnames (chainsPTIRT[[1]]))],
                chainsPTIRT[[2]][,grep("X",  colnames (chainsPTIRT[[2]]))])

Deviance <- c (chainsPTIRT[[1]][,grep("deviance",  colnames (chainsPTIRT[[1]]))],
               chainsPTIRT[[2]][,grep("deviance",  colnames (chainsPTIRT[[2]]))])


item.labels <- c("Party Sponsor",   # substanceSponsorParty
                 "Party Leader's Vote Position",    #partyLeaderPositionInMedia
                 "Liberated Vote",   # liberate
                 "Conjectured Party Position",   # PartyPositionSuggestedByMedia
                 "Policy Elaboration",  # PolicyElaboration
                 "Procedural Vote",   # proceduralVoteLeaderRequest
                 "Obstruction",  #  obstruction
                 "Leadership Floor Speech",  # partyLeaderIntroducesPosition
                 "Divided Speech",   # dividedSpeech
                 "Party President's Position",  #nationalPartyPresidentPositionInMedia
                 "Party Leader's Bill Position", # partyLeaderBillPositionInMedia
                 "Leadership Cohesiveness") # LeaderStrongVotePositionTaking



# Party pressure scores
X.pt    <- X.e[,grep(",1]", colnames(X.e))]
X.pmdb  <- X.e[,grep(",2]", colnames(X.e))]
X.psdb  <- X.e[,grep(",3]", colnames(X.e))]
X.dem   <- X.e[,grep(",4]", colnames(X.e))]
X.pp    <- X.e[,grep(",5]", colnames(X.e))]
X.pr    <- X.e[,grep(",6]", colnames(X.e))]

# Graphs with uncertainty about pressure scores
# Figure C2 Appendix (PT)

# pdf (paste0(graphPath, "pressurePT.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pt)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PT")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.pt), colMeans(X.pt))
        , pch=20)
segments (x0=1:ncol(X.pt), x1=1:ncol(X.pt)
          , y0=apply (X.pt, 2, quantile, 0.1)
          , y1=apply (X.pt, 2, quantile, 0.9))
# dev.off()

# Figure C2 Appendix (PMDB)
# pdf (paste0(graphPath, "pressurePMDB.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pmdb)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PMDB")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.pmdb), colMeans(X.pmdb))
        , pch=20)
segments (x0=1:ncol(X.pmdb), x1=1:ncol(X.pmdb)
          , y0=apply (X.pmdb, 2, quantile, 0.1)
          , y1=apply (X.pmdb, 2, quantile, 0.9))
# dev.off()

# Figure C2 Appendix (PSDB)
# pdf (paste0(graphPath, "pressurePSDB.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.psdb)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PSDB")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.psdb), colMeans(X.psdb))
        , pch=20)
segments (x0=1:ncol(X.psdb), x1=1:ncol(X.psdb)
          , y0=apply (X.psdb, 2, quantile, 0.1)
          , y1=apply (X.psdb, 2, quantile, 0.9))
# dev.off()

# Figure C2 Appendix (DEM)
# pdf (paste0(graphPath, "pressureDEM.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.dem)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="DEM")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.dem), colMeans(X.dem))
        , pch=20)
segments (x0=1:ncol(X.dem), x1=1:ncol(X.dem)
          , y0=apply (X.dem, 2, quantile, 0.1)
          , y1=apply (X.dem, 2, quantile, 0.9))
# dev.off()

# Figure C2 Appendix (PP)
# pdf (paste0(graphPath, "pressurePP.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pp)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PP")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.pp), colMeans(X.pp))
        , pch=20)
segments (x0=1:ncol(X.pp), x1=1:ncol(X.pp)
          , y0=apply (X.pp, 2, quantile, 0.1)
          , y1=apply (X.pp, 2, quantile, 0.9))
# dev.off()

# Figure C2 Appendix (PR)
# pdf (paste0(graphPath, "pressurePR.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pr)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PR")
axis (2)
axis (1)
points (xy.coords (1:ncol(X.pr), colMeans(X.pr))
        , pch=20)
segments (x0=1:ncol(X.pr), x1=1:ncol(X.pr)
          , y0=apply (X.pr, 2, quantile, 0.1)
          , y1=apply (X.pr, 2, quantile, 0.9))
# dev.off()

# Pressure discrimination parameters
Beta.pt    <- Beta.e[,grep(",1]", colnames(Beta.e))]
Beta.pmdb  <- Beta.e[,grep(",2]", colnames(Beta.e))]
Beta.psdb  <- Beta.e[,grep(",3]", colnames(Beta.e))]
Beta.dem  <- Beta.e[,grep(",4]", colnames(Beta.e))]
Beta.pp   <- Beta.e[,grep(",5]", colnames(Beta.e))]
Beta.pr    <- Beta.e[,grep(",6]", colnames(Beta.e))]

# Pressure difficulty parameters
Alpha.pt    <- Alpha.e[,grep(",1]", colnames(Alpha.e))]
Alpha.pmdb  <- Alpha.e[,grep(",2]", colnames(Alpha.e))]
Alpha.psdb  <- Alpha.e[,grep(",3]", colnames(Alpha.e))]
Alpha.dem   <- Alpha.e[,grep(",4]", colnames(Alpha.e))]
Alpha.pp    <- Alpha.e[,grep(",5]", colnames(Alpha.e))]
Alpha.pr    <- Alpha.e[,grep(",6]", colnames(Alpha.e))]



# Difficulty parameters
# Figure C1.b Appendix 
# pdf (paste0(graphPath, "fig8.pdf"), h=10, w=10)
par (las=2, mar=c(12,3,1,1), cex.axis=1)
plot (colMeans (Alpha.pt), ylim=c(0,40), pch=19, col="red"
      , axes=F, xlab="", ylab="Posterior Alpha")
segments (x0=c(1:ncol(Alpha.pt)), x1=c(1:ncol(Alpha.pt))
          , y0=colMeans (Alpha.pt) + apply (Alpha.pt, 2, sd)
          , y1=colMeans (Alpha.pt) - apply (Alpha.pt, 2, sd)
          , col="red")

points (xy.coords (c(1:12) + 0.1, colMeans (Alpha.pmdb)), col="green", pch=19)
segments (x0=c(1:ncol(Alpha.pmdb)) + 0.1, x1=c(1:ncol(Alpha.pmdb) + 0.1)
          , y0=colMeans (Alpha.pmdb) + apply (Alpha.pmdb, 2, sd)
          , y1=colMeans (Alpha.pmdb) - apply (Alpha.pmdb, 2, sd)
          , col="green")

points (xy.coords (c(1:12) - 0.1, colMeans (Alpha.psdb)), col="purple", pch=19)
segments (x0=c(1:ncol(Alpha.psdb)) - 0.1, x1=c(1:ncol(Alpha.psdb) - 0.1)
          , y0=colMeans (Alpha.psdb) + apply (Alpha.psdb, 2, sd)
          , y1=colMeans (Alpha.psdb) - apply (Alpha.psdb, 2, sd)
          , col="purple")

points (xy.coords (c(1:12) - 0.2, colMeans (Alpha.dem)), col="black", pch=19)
segments (x0=c(1:ncol(Alpha.dem)) - 0.2, x1=c(1:ncol(Alpha.dem) - 0.2)
          , y0=colMeans (Alpha.dem) + apply (Alpha.dem, 2, sd)
          , y1=colMeans (Alpha.dem) - apply (Alpha.dem, 2, sd)
          , col="black")

points (xy.coords (c(1:12) + 0.2, colMeans (Alpha.pp)), col="orange", pch=19)
segments (x0=c(1:ncol(Alpha.pp)) + 0.2, x1=c(1:ncol(Alpha.pp) + 0.2)
          , y0=colMeans (Alpha.pp) + apply (Alpha.pp, 2, sd)
          , y1=colMeans (Alpha.pp) - apply (Alpha.pp, 2, sd)
          , col="orange")

points (xy.coords (c(1:12) - 0.3, colMeans (Alpha.pr)), col="gray", pch=19)
segments (x0=c(1:ncol(Alpha.pr)) - 0.3, x1=c(1:ncol(Alpha.pr) - 0.3)
          , y0=colMeans (Alpha.pr) + apply (Alpha.pr, 2, sd)
          , y1=colMeans (Alpha.pr) - apply (Alpha.pr, 2, sd)
          , col="gray")

abline (h=0, lty=3)
abline (v=c(1:12)+0.45, lty=3)
axis (2)
axis (1, at=c(1:ncol(Alpha.pt)),
      labels=item.labels,
      cex=1)
legend ("top",
        legend=c("PT","PMDB","PSDB","DEM","PP","PR"),
        pch=19,
        col=c("red","green","purple","black","orange","grey"))
# dev.off()



# Discrimination parameters
# Figure C1.a Appendix 
# pdf (paste0(graphPath, "fig7.pdf"), h=10, w=10)
par (las=2, mar=c(12,3,1,1), cex.axis=1)
plot (colMeans (Beta.pt), ylim=c(-5,6), pch=19, col="red"
      , axes=F, xlab="", ylab="Posterior beta")
segments (x0=c(1:ncol(Beta.pt)), x1=c(1:ncol(Beta.pt))
          , y0=colMeans (Beta.pt) + apply (Beta.pt, 2, sd)
          , y1=colMeans (Beta.pt) - apply (Beta.pt, 2, sd)
          , col="red")

points (xy.coords (c(1:12) + 0.1, colMeans (Beta.pmdb)), col="green", pch=19)
segments (x0=c(1:ncol(Beta.pmdb)) + 0.1, x1=c(1:ncol(Beta.pmdb) + 0.1)
          , y0=colMeans (Beta.pmdb) + apply (Beta.pmdb, 2, sd)
          , y1=colMeans (Beta.pmdb) - apply (Beta.pmdb, 2, sd)
          , col="green")

points (xy.coords (c(1:12) - 0.1, colMeans (Beta.psdb)), col="purple", pch=19)
segments (x0=c(1:ncol(Beta.psdb)) - 0.1, x1=c(1:ncol(Beta.psdb) - 0.1)
          , y0=colMeans (Beta.psdb) + apply (Beta.psdb, 2, sd)
          , y1=colMeans (Beta.psdb) - apply (Beta.psdb, 2, sd)
          , col="purple")

points (xy.coords (c(1:12) - 0.2, colMeans (Beta.dem)), col="black", pch=19)
segments (x0=c(1:ncol(Beta.dem)) - 0.2, x1=c(1:ncol(Beta.dem) - 0.2)
          , y0=colMeans (Beta.dem) + apply (Beta.dem, 2, sd)
          , y1=colMeans (Beta.dem) - apply (Beta.dem, 2, sd)
          , col="black")

points (xy.coords (c(1:12) + 0.2, colMeans (Beta.pp)), col="orange", pch=19)
segments (x0=c(1:ncol(Beta.pp)) + 0.2, x1=c(1:ncol(Beta.pp) + 0.2)
          , y0=colMeans (Beta.pp) + apply (Beta.pp, 2, sd)
          , y1=colMeans (Beta.pp) - apply (Beta.pp, 2, sd)
          , col="orange")

points (xy.coords (c(1:12) - 0.3, colMeans (Beta.pr)), col="grey", pch=19)
segments (x0=c(1:ncol(Beta.pr)) - 0.3, x1=c(1:ncol(Beta.pr) - 0.3)
          , y0=colMeans (Beta.pr) + apply (Beta.pr, 2, sd)
          , y1=colMeans (Beta.pr) - apply (Beta.pr, 2, sd)
          , col="grey")


abline (h=0, lty=3)
abline (v=c(1:12)+0.3, lty=3)
axis (2)
axis (1, at=c(1:ncol(Beta.pt)),
      labels=item.labels,
      cex=1)
legend ("bottomright", 
        legend=c("PT","PMDB","PSDB","DEM","PP","PR"),
        pch=19,
        col=c("red","green","purple","black","orange","grey"))
# dev.off ()


## Figure C3a (Appendix) - top row
# pdf (paste0(graphPath, "fig9.pdf"), h=6, w=16)  
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2, las=0)
hist (colMeans (X.pt)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PT", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.pmdb)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PMDB", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.pp)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PP", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
# dev.off()

## Figure C3a (Appendix) - bottom row
# pdf (paste0(graphPath, "fig11.pdf"), h=6, w=16)
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2, las=0)
hist (colMeans (X.pr)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PR", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.psdb)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PSDB", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.dem)
      , ylim=c(0,140)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="DEM", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
# dev.off()



#### Plotting directional party pressure ####

# Directional party pressure graphs for Brazil
load("brazil_rollcall_long_format_mean.Rdata")

party <- unlist (strsplit (pressureVector$id, split="-"))[c(TRUE,FALSE)]

## Figure C3b (Appendix) - top row
# pdf (paste0 (graphPath, "fig10.pdf"), h=6, w=16)  
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2)
hist (pressureVector$pressure[party=="PT"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=25
      , main="PT", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressureVector$pressure[party=="PMDB"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="PMDB", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressureVector$pressure[party=="PP"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="PP", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
# dev.off()

## Figure C3b (Appendix) - bottom row
# pdf (paste0 (graphPath, "fig12.pdf"), h=6, w=16)  
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2)
hist (pressureVector$pressure[party=="PR"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=25
      , main="PR", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressureVector$pressure[party=="PSDB"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="PSDB", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressureVector$pressure[party=="DEM"]
      , ylim=c(0,100)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="DEM", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
# dev.off()




###############################
#### Mexico party pressure ####
###############################
rm (list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()  

load("mexico_ptirt_mcmc.Rda")

# Check convergence
gD <- gelman.diag (chainsPTIRT, multivariate=F) 
summary (gD[[1]][,1]) 

Beta.e   <- rbind (chainsPTIRT[[1]][,grep("beta",  colnames (chainsPTIRT[[1]]))],
                   chainsPTIRT[[2]][,grep("beta",  colnames (chainsPTIRT[[2]]))])

Alpha.e   <- rbind (chainsPTIRT[[1]][,grep("alpha",  colnames (chainsPTIRT[[1]]))],
                    chainsPTIRT[[2]][,grep("alpha",  colnames (chainsPTIRT[[2]]))])

X.e   <- rbind (chainsPTIRT[[1]][,grep("X",  colnames (chainsPTIRT[[1]]))],
                chainsPTIRT[[2]][,grep("X",  colnames (chainsPTIRT[[2]]))])

Deviance <- c (chainsPTIRT[[1]][,grep("deviance",  colnames (chainsPTIRT[[1]]))],
               chainsPTIRT[[2]][,grep("deviance",  colnames (chainsPTIRT[[2]]))])

X.pri  <- X.e[,grep(",1]", colnames(X.e))]
X.pan  <- X.e[,grep(",2]", colnames(X.e))]
X.prd  <- X.e[,grep(",3]", colnames(X.e))]

Beta.pri <- Beta.e[,grep(",1]", colnames (Beta.e))]
Beta.pan <- Beta.e[,grep(",2]", colnames (Beta.e))]
Beta.prd <- Beta.e[,grep(",3]", colnames (Beta.e))]

Alpha.pri <- Alpha.e[,grep(",1]", colnames (Alpha.e))]
Alpha.pan <- Alpha.e[,grep(",2]", colnames (Alpha.e))]
Alpha.prd <- Alpha.e[,grep(",3]", colnames (Alpha.e))]

# Graphs with uncertainty on party pressure estimates
## Figure C4 (Appendix) - PRI
par (mfrow=c(1,1))
# pdf (paste0(graphPath, "pressurePRI.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pri)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PRI")
axis (2, at=c(0,1,2,3), labels=c(0,1,2,3))
axis (1, at=seq(0,270,10), labels=seq(0,270,10))
points (xy.coords (1:ncol(X.pri), colMeans(X.pri))
        , pch=20)
segments (x0=1:ncol(X.pri), x1=1:ncol(X.pri)
          , y0=apply (X.pri, 2, quantile, 0.1)
          , y1=apply (X.pri, 2, quantile, 0.9))
# dev.off()

## Figure C4 (Appendix) - PAN
# pdf (paste0(graphPath, "pressurePAN.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.pan)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PAN")
axis (2, at=c(0,1,2,3), labels=c(0,1,2,3))
axis (1, at=seq(0,270,10), labels=seq(0,270,10))
points (xy.coords (1:ncol(X.pan), colMeans(X.pan))
        , pch=20)
segments (x0=1:ncol(X.pan), x1=1:ncol(X.pan)
          , y0=apply (X.pan, 2, quantile, 0.1)
          , y1=apply (X.pan, 2, quantile, 0.9))
# dev.off()

## Figure C4 (Appendix) - PRD
# pdf (paste0(graphPath, "pressurePRD.pdf"), h=6, w=8)
par (las=2, mar=c(4,4,3,1), cex.main=2)
plot (x=c(1,ncol(X.prd)), y=c(0,3), type="n"
      , axes=F
      , xlab="Vote"
      , ylab="Estimated party pressure"
      , main="PRD")
axis (2, at=c(0,1,2,3), labels=c(0,1,2,3))
axis (1, at=seq(0,270,10), labels=seq(0,270,10))
points (xy.coords (1:ncol(X.prd), colMeans(X.prd))
        , pch=20)
segments (x0=1:ncol(X.prd), x1=1:ncol(X.prd)
          , y0=apply (X.prd, 2, quantile, 0.1)
          , y1=apply (X.prd, 2, quantile, 0.9))
# dev.off()




# Figure C6a (Appendix)
# pdf (paste0 (graphPath, "fig13.pdf"), h=6, w=16)  
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2)
hist (colMeans (X.pri), ylim=c(0,70)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PRI", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.pan), ylim=c(0,70)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PAN", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
hist (colMeans (X.prd), ylim=c(0,70)
      , xlim=c(0.5,2.2)
      , breaks=15
      , main="PRD", xlab="")
mtext (1, line=2.5, text="Pressure scores", cex=1)
# dev.off()




item.labels <- c("Party President's Position", # "PRI.presidentPartyPosition"
                 "Party Leader's Bill Position", # "PRI.partyLeaderPosition"
                 "Committee Bill", # "PRI.committeeBill"
                 "Party Leader's Vote Position", # "PRI.PartyLeaderVotePosition",
                 "Leadership Cohesiveness", # "PRI.LeaderVotePositionTaking"
                 "Leadership Speech", #"PRI.leadershipSpeech"
                 "Delegate Threat", # "PRI.delegateThreat"
                 "Delegate Position", # "PRI.delegatePosition"
                 "Position Declaration", # "PRI.positionDeclaration"
                 "Weak Leadership Position") # PRI.LeaderWeakPosition

# Beta plots
# Figure C5a (Appendix)
# pdf (paste0(graphPath, "fig15.pdf"), h=10, w=10)
par (las=2, mar=c(12,3,1,1), cex.axis=1, mfrow=c(1,1))
plot ((colMeans (Beta.pri)), ylim=c(-2.5,8), pch=19, col="green"
      , axes=F, xlab="", ylab="Posterior beta", type="n")
points (xy.coords (c(1:10) - 0.2, colMeans (Beta.pri)), col="green", pch=19)
segments (x0=c(1:ncol(Beta.pri)) - 0.2, x1=c(1:ncol(Beta.pri) - 0.2)
          , y0=colMeans (Beta.pri) + apply (Beta.pri, 2, sd)
          , y1=colMeans (Beta.pri) - apply (Beta.pri, 2, sd)
          , col="green")
points (xy.coords (c(1:10), colMeans (Beta.pan)), col="blue", pch=19)
segments (x0=c(1:ncol(Beta.pan)), x1=c(1:ncol(Beta.pan))
          , y0=colMeans (Beta.pan) + apply (Beta.pan, 2, sd)
          , y1=colMeans (Beta.pan) - apply (Beta.pan, 2, sd)
          , col="blue")
points (xy.coords (c(1:10) + 0.2, colMeans (Beta.prd)), col="red", pch=19)
segments (x0=c(1:ncol(Beta.prd)) + 0.2, x1=c(1:ncol(Beta.prd) + 0.2)
          , y0=colMeans (Beta.prd) + apply (Beta.prd, 2, sd)
          , y1=colMeans (Beta.prd) - apply (Beta.prd, 2, sd)
          , col="red")
abline (h=0, lty=3)
axis (2)
axis (1, at=c(1:ncol(Beta.pri)),
      labels=item.labels,
      cex=1)
legend ("topright", 
        legend=c("PRI","PAN","PRD"),
        pch=19,
        col=c("green","blue","red"))
# dev.off()

# Alpha plots
# Figure C5b (Appendix)
# pdf (paste0(graphPath, "fig16.pdf"), h=10, w=10)
par (las=2, mar=c(12,3,1,1), cex.axis=1)
plot (colMeans (Alpha.pri), ylim=c(0,30), pch=19, col="green"
      , axes=F, xlab="", ylab="Posterior Alpha", type="n")
points (xy.coords (c(1:10) - 0.2, colMeans (Alpha.pri))
        , col="green", pch=19)
segments (x0=c(1:ncol(Alpha.pri)) - 0.2, x1=c(1:ncol(Alpha.pri)) - 0.2
          , y0=colMeans (Alpha.pri) + apply (Alpha.pri, 2, sd)
          , y1=colMeans (Alpha.pri) - apply (Alpha.pri, 2, sd)
          , col="green")
points (xy.coords (c(1:10), colMeans (Alpha.pan))
        , col="blue", pch=19)
segments (x0=c(1:ncol(Alpha.pan)), x1=c(1:ncol(Alpha.pan))
          , y0=colMeans (Alpha.pan) + apply (Alpha.pan, 2, sd)
          , y1=colMeans (Alpha.pan) - apply (Alpha.pan, 2, sd)
          , col="blue")
points (xy.coords (c(1:10) + 0.2, colMeans (Alpha.prd)), col="red", pch=19)
segments (x0=c(1:ncol(Alpha.prd)) + 0.2, x1=c(1:ncol(Alpha.prd) + 0.2)
          , y0=colMeans (Alpha.prd) + apply (Alpha.prd, 2, sd)
          , y1=colMeans (Alpha.prd) - apply (Alpha.prd, 2, sd)
          , col="red")
abline (h=0, lty=3)
axis (2)
axis (1, at=c(1:ncol(Beta.pri)),
      labels=item.labels,
      cex=1)
# dev.off()


## Directional party pressure
PRI.pressure.index <- apply (X.pri, 2, mean)
PAN.pressure.index <- apply (X.pan, 2, mean)
PRD.pressure.index <- apply (X.prd, 2, mean)

# Load PRI, PAN, PRD leader positions 
load (file="mexico_leader_decisions.Rda")
PRILeaderDecision <- objects2keep[[1]]
PANLeaderDecision <- objects2keep[[2]]
PRDLeaderDecision <- objects2keep[[3]]
rm (objects2keep)

pressureMatrix <- rbind (PRI.pressure.index, PAN.pressure.index, PRD.pressure.index)
pressureMatrix.exp <-  rbind (pressureMatrix, rep(0,ncol(pressureMatrix)), rep(0,ncol(pressureMatrix)), rep(0,ncol(pressureMatrix)), rep(0,ncol(pressureMatrix)), rep(0,ncol(pressureMatrix))) # pressureMatrix # 
positionMatrix <- rbind(PRILeaderDecision    # PRI.hi.std.LeaderDecision contains more "unknown" party positions than PRILeaderDecision
                        , PANLeaderDecision  # same as above
                        , PRDLeaderDecision)
positionMatrix.exp <- rbind (positionMatrix, rep(0,ncol(positionMatrix)), rep(0,ncol(positionMatrix)), rep(0,ncol(positionMatrix)), rep(0,ncol(positionMatrix)), rep(0,ncol(positionMatrix))) # positionMatrix #

pressure.for <- pressureMatrix.exp*positionMatrix.exp

# Figure C6b (Appendix)
# pdf (paste0(graphPath, "fig14.pdf"), h=6, w=16) 
par (mfrow=c(1,3), mar=c(4,3,2,1), cex.main=2, las=0)
hist (pressure.for[1,]
      , ylim=c(0,90)
      , xlim=c(-2.2,2.2)
      , breaks=25
      , main="PRI", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressure.for[2,]
      , ylim=c(0,90)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="PAN", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
hist (pressure.for[3,]
      , ylim=c(0,90)
      , xlim=c(-2.2,2.2)
      , breaks=15
      , main="PRD", xlab="")
mtext (1, line=2.5, text="Directional pressure", cex=1)
# dev.off()  # 

par (mfrow=c(1,1))

################################
#### Report model summaries ####
################################


files2load <- c("objects_2MPE_Mexico.Rdata",
                "objects_1MPE_Mexico.Rdata")

variablesNames <- list ( TwoMPE=c()
                       , OneMPE=c())

niceNames <- c("PRI pressure","PAN pressure","PRD pressure","Pressure x PR", "PR"
               , "Pressure x Poverty", "Poverty"
               , "Pressure x Year","Year"
               , "Pressure x Victory Margin", "Victory Margin"
               , "Pressure x State Legislators", "State Legislators"
               , "Pressure x Cabinet Posts", "Cabinet Posts"
               , "Pressure x Presidential Position", "Presidential Position"
               , "Pressure x Committee Chair",  "Committee Chair"
               , "Pressure x PCB Request","PCB Request")


## Mexico table, with standard deviations
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  load (files2load[i])
  
  # Theta  <- rbind (chainsComp[[1]][, grep("theta",  colnames (chainsComp[[1]]))]
  #                  , chainsComp[[2]][, grep("theta",  colnames (chainsComp[[2]]))])
  # Delta  <- rbind (chainsComp[[1]][, grep("delta",  colnames (chainsComp[[1]]))]
  #                  , chainsComp[[2]][, grep("delta",  colnames (chainsComp[[2]]))])
  Lambda <- items2export$lambda
  mean.lambda <- colMeans (Lambda)
  sd.lambda   <- apply (Lambda, 2, sd)
  star <- ifelse (abs (mean.lambda/sd.lambda) > 1.96, "*", "")
  mean.lambda <- paste0 (round (mean.lambda, 2), star)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, round (sd.lambda, 2))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}
rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames



## Mexico table, with 99% credible intervals
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  load (files2load[i])
  Lambda <- items2export$lambda
  mean.lambda <- colMeans (Lambda)
  sd.lambda   <- apply (Lambda, 2, sd)
  lo.quant   <- apply (Lambda, 2, quantile, prob=0.005)
  hi.quant   <- apply (Lambda, 2, quantile, prob=0.995)
  mean.lambda <- round (mean.lambda, 2)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, paste0("[", round (lo.quant, 2), ":", round (hi.quant, 2), "]"))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}
rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames

rearranged.all.means.lambda <- all.means.lambda[c(1,2,3,5,7,11,17,21,9,13,15,19,4,6,10,16,20,8,12,14,18), ]

# Table H2 (Appendix) 
# latex ref t:mexicoTable2-appendix
xtable (rearranged.all.means.lambda
        , caption="The effects of specific characteristics on the propensity to vote Aye ($\\hat{\\mu}$) in Brazil and Mexico (complementary results for Tables 5 and 6)."
        , label="T:mexicoTable2")

# The following quantities of interest need to be added by hand to t:mexicoTable2
load (files2load[1])
probVoteDisagree  <- items2export$voteDis
probVoteAgree     <- items2export$voteAgr
probVoteUncertain <- items2export$voteUnc
  
round (quantile (c(probVoteAgree, probVoteAgree), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteUncertain, probVoteUncertain), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteDisagree, probVoteDisagree), prob=c(0.005, 0.5, 0.995)), 3)
rm (items2export)



# Table H4 (Appendix) 
# TABLE t:mexicoTable2_appendix_posterior_mean,
# based on one draw, the mean of the posterior distribution of party pressure
# Load 1mpe files
load ("mexico_pressure_final_result_sample_meana.Rda")
process.1 <- chainsPress; rm(chainsPress)

load ("mexico_pressure_final_result_sample_meanb.Rda")
process.2 <- chainsPress; rm(chainsPress)

allChains <- as.mcmc.list (list (process.1, process.2))
gelman.Comp <- gelman.diag (allChains, multivariate=F) # 
print (max (gelman.Comp[[1]][,1], na.rm = TRUE))

lambda <- deviance <- c()
lm1 <- process.1[,grep("lambda", colnames(process.1))]
dev1 <- process.1[,grep("deviance", colnames (process.1))]
lm2 <- process.2[,grep("lambda", colnames(process.2))]
dev2 <- process.2[,grep("deviance", colnames (process.2))]
  
lm <- rbind (lm1, lm2)
dev  <- c(dev1, dev2)
  
lambda.1mpe <- rbind (lambda, lm)
deviance.1mpe <- c(deviance, dev)

# Load 2mpe files
rm (process.1, process.2)
load ("mexico_competing_final_result_sample_meana.Rda")
process.1 <- chainsComp; rm(chainsComp)

load ("mexico_competing_final_result_sample_meanb.Rda")
process.2 <- chainsComp; rm(chainsComp)

allChains <- as.mcmc.list (list (process.1, process.2))
gelman.Comp <- gelman.diag (allChains, multivariate=F) # 
print (max (gelman.Comp[[1]][,1], na.rm = TRUE))

lambda <- deviance <- c()
lm1 <- process.1[,grep("lambda", colnames(process.1))]
dev1 <- process.1[,grep("deviance", colnames (process.1))]
agr1 <- process.1[,grep("probVoteAgree", colnames(process.1))]
dis1 <- process.1[,grep("probVoteDisagree", colnames(process.1))]
unc1 <- process.1[,grep("probVoteUncertain", colnames(process.1))]
lm2 <- process.2[,grep("lambda", colnames(process.2))]
dev2 <- process.2[,grep("deviance", colnames (process.2))]
agr2 <- process.2[,grep("probVoteAgree", colnames(process.2))]
dis2 <- process.2[,grep("probVoteDisagree", colnames(process.2))]
unc2 <- process.2[,grep("probVoteUncertain", colnames(process.2))]

lm <- rbind (lm1, lm2)
dev  <- c(dev1, dev2)
agr  <- c(agr1, agr2)
dis  <- c(dis1, dis2) 
unc  <- c(unc1, unc2)

lambda.2mpe <- rbind (lambda, lm)
deviance.2mpe <- c(deviance, dev)
probVoteDisagree  <- dis
probVoteAgree     <- agr
probVoteUncertain <- unc

Lambda <- list (lambda.2mpe, lambda.1mpe)
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  mean.lambda <- colMeans (Lambda[[i]])
  sd.lambda   <- apply (Lambda[[i]], 2, sd)
  lo.quant   <- apply (Lambda[[i]], 2, quantile, prob=0.005)
  hi.quant   <- apply (Lambda[[i]], 2, quantile, prob=0.995)
  # star <- ifelse (abs (mean.lambda/sd.lambda) > 1.96, "*", "")
  mean.lambda <- round (mean.lambda, 2)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, paste0("[", round (lo.quant, 2), ":", round (hi.quant, 2), "]"))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}
rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames

rearranged.all.means.lambda <- all.means.lambda[c(1,2,3,5,7,11,17,21,9,13,15,19,4,6,10,16,20,8,12,14,18), ]


# Table H4 (Appendix) 
xtable (rearranged.all.means.lambda
        , caption="The effects of specific characteristics on the propensity to vote Aye ($\\hat{\\mu}$) in Mexico (results based on party pressure point estimate, which is the mean of the posterior distribution of the PTIRT model)."
        , label="T:mexicoTable_appendix")

# The following quantities of interest need to be added by hand to t:mexicoTable2
round (quantile (c(probVoteAgree, probVoteAgree), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteUncertain, probVoteUncertain), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteDisagree, probVoteDisagree), prob=c(0.005, 0.5, 0.995)), 3)







############
###Brazil###
############
rm (list=ls())


# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()   #"/Users/grosas/Dropbox/PartyEffects/MacPro/JOP RnR/newGraphs/"

files2load <- c("objects_2MPE_Brazil.Rdata",
                "objects_1MPE_Brazil.Rdata")

variablesNames <- list ( TwoMPE=c(),
                         OneMPE=c())

niceNames <- c("PT directional pressure","PSDB directional pressure"
               , "DEM directional pressure","PMDB directional pressure"
               , "PR directional pressure","PP directional pressure"
               , "Directional pressure x Seniority","Seniority"
               , "Directional pressure x District income","District income"
               , "Directional pressure x Commmittee chair", "Committee chair"
               , "Directional pressure x Percentile rank", "Percentile rank"
               , "Directional pressure x District magnitude","District magnitude"
               , "Directional pressure x Year","Year"
               , "Directional pressure x Procedural vote","Procedural vote"
               , "Directional pressure x Cabinet post", "Cabinet post"
               , "Directional pressure x Prez position", "Prez position"
               , "Directional pressure x Party agreements","Party agreements")

## Brazil table, with standard deviations
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  load (files2load[i])
  Lambda <- items2export$lambda
  mean.lambda <- colMeans (Lambda)
  sd.lambda   <- apply (Lambda, 2, sd)
  star <- ifelse (abs (mean.lambda/sd.lambda) > 1.96, "*", "")
  mean.lambda <- paste0 (round (mean.lambda, 2), star)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, round (sd.lambda, 2))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}

rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames


## Brazil table, with 99% credible intervals
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  load (files2load[i])
  Lambda <- items2export$lambda
  mean.lambda <- colMeans (Lambda)
  sd.lambda   <- apply (Lambda, 2, sd)
  mean.lambda <- round (mean.lambda, 2)
  lo.quant   <- apply (Lambda, 2, quantile, prob=0.005)
  hi.quant   <- apply (Lambda, 2, quantile, prob=0.995)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, paste0("[", round (lo.quant, 2), ":", round (hi.quant, 2), "]"))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}

rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames

rearranged.all.means.lambda.BZ <- all.means.lambda[c(1,2,3,4,5,6,10,14,16,20,24,26,8,12,18,22,9,13,15,19,23,25,7,11,17,21),]

# Table H1 (Appendix)
# latex ref t:competingBrazil-appendix
xtable (rearranged.all.means.lambda.BZ
        , caption="Party effects in Brazil"
        , label="t:competingBrazil")

# The following quantities of interest need to be added by hand to t:competingBrazil
load (files2load[1])
probVoteDisagree  <- items2export$voteDis
probVoteAgree     <- items2export$voteAgr
probVoteUncertain <- items2export$voteUnc

round (quantile (c(probVoteAgree, probVoteAgree), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteUncertain, probVoteUncertain), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteDisagree, probVoteDisagree), prob=c(0.005, 0.5, 0.995)), 3)




# latex ref t:competingBrazil_Appendix_posterior_mean
# based on one draw, the mean of the posterior distribution of party pressure
# Load 1mpe files
rm (process.1, process.2)
load ("brazil_pressure_final_result_mean_a.Rda")
process.1 <- as.mcmc(chainsPress); rm(chainsPress)

load ("brazil_pressure_final_result_mean_b.Rda")
process.2 <- as.mcmc(chainsPress); rm(chainsPress)

allChains <- list (process.1, process.2)
gelman.Comp <- gelman.diag (allChains, multivariate=F) # 
print (max (gelman.Comp[[1]][,1], na.rm = TRUE))

lambda <- deviance <- c()
lm1 <- process.1[,grep("lambda", colnames(process.1))]
dev1 <- process.1[,grep("deviance", colnames (process.1))]
lm2 <- process.2[,grep("lambda", colnames(process.2))]
dev2 <- process.2[,grep("deviance", colnames (process.2))]

lm <- rbind (lm1, lm2)
dev  <- c(dev1, dev2)

lambda.1mpe <- rbind (lambda, lm)
deviance.1mpe <- c(deviance, dev)


# Load 2mpe files
rm (process.1, process.2)
load ("brazil_competing_final_result_sample_meana.Rda")
process.1 <- as.mcmc(chainsComp); rm(chainsComp)

load ("brazil_competing_final_result_sample_meanb.Rda")
process.2 <- as.mcmc(chainsComp); rm(chainsComp)

allChains <- as.mcmc.list (list (process.1, process.2))
gelman.Comp <- gelman.diag (allChains, multivariate=F) # 
print (max (gelman.Comp[[1]][,1], na.rm = TRUE))

lambda <- deviance <- c()
lm1 <- process.1[,grep("lambda", colnames(process.1))]
dev1 <- process.1[,grep("deviance", colnames (process.1))]
agr1 <- process.1[,grep("probVoteAgree", colnames(process.1))]
dis1 <- process.1[,grep("probVoteDisagree", colnames(process.1))]
unc1 <- process.1[,grep("probVoteUncertain", colnames(process.1))]
lm2 <- process.2[,grep("lambda", colnames(process.2))]
dev2 <- process.2[,grep("deviance", colnames (process.2))]
agr2 <- process.2[,grep("probVoteAgree", colnames(process.2))]
dis2 <- process.2[,grep("probVoteDisagree", colnames(process.2))]
unc2 <- process.2[,grep("probVoteUncertain", colnames(process.2))]

lm <- rbind (lm1, lm2)
dev  <- c(dev1, dev2)
agr  <- c(agr1, agr2)
dis  <- c(dis1, dis2) 
unc  <- c(unc1, unc2)

lambda.2mpe <- rbind (lambda, lm)
deviance.2mpe <- c(deviance, dev)
probVoteDisagree  <- dis
probVoteAgree     <- agr
probVoteUncertain <- unc


# Arrange data into table
Lambda <- list (lambda.2mpe, lambda.1mpe)
all.means.lambda <- c()
columnNames <- c()
for (i in 1:2){
  mean.lambda <- colMeans (Lambda[[i]])
  sd.lambda   <- apply (Lambda[[i]], 2, sd)
  lo.quant   <- apply (Lambda[[i]], 2, quantile, prob=0.005)
  hi.quant   <- apply (Lambda[[i]], 2, quantile, prob=0.995)
  # star <- ifelse (abs (mean.lambda/sd.lambda) > 1.96, "*", "")
  mean.lambda <- round (mean.lambda, 2)
  all.means.lambda <- cbind (all.means.lambda, mean.lambda, paste0("[", round (lo.quant, 2), ":", round (hi.quant, 2), "]"))
  columnNames <- c(columnNames, rep (names (variablesNames)[i], 2))
}
rownames (all.means.lambda) <- niceNames
colnames (all.means.lambda) <- columnNames

rearranged.all.means.lambda.BZ <- all.means.lambda[c(1,2,3,4,5,6,10,14,16,20,24,26,8,12,18,22,9,13,15,19,23,25,7,11,17,21),]

# Table H3 (Appendix)
xtable (rearranged.all.means.lambda.BZ
        , caption="Party effects in Brazil, results based on party pressure point estimates, omitting estimation uncertainty from the party pressure measurement model."
        , label="t:competingBrazil_Appendix_posterior_mean")


# The following quantities of interest need to be added by hand to t:competingBrazil_Appendix_posterior_mean
round (quantile (c(probVoteAgree, probVoteAgree), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteUncertain, probVoteUncertain), prob=c(0.005, 0.5, 0.995)), 3)
round (quantile (c(probVoteDisagree, probVoteDisagree), prob=c(0.005, 0.5, 0.995)), 3)



#####################################
###plot distribution of ideal point##
#####################################
##############
### Brazil ###
##############
rm (list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()  
load (file="brazil_rollcall_long_format_mean.Rdata")  

new.party.color <- as.numeric (unlist( by (molten.party.member$party.no, molten.party.member$final.leg.no, unique)))
new.party.color <- car::recode(new.party.color, " 1='red'
                               ; 2='purple'
                               ; 3='orange'
                               ; 4='yellow'
                               ; 5='black'
                               ; 6='blue'
                               ; 7='gray'")


##Plot CJR 
load (file="brazil_CJR_auxiliary_output.Rda")


chains <- chainsCJR

th1 <- chains[[1]][,grep("theta", colnames(chains[[1]]))]
dt1 <- chains[[1]][,grep("delta", colnames(chains[[1]]))]
bt1 <- chains[[1]][,grep("beta", colnames(chains[[1]]))]
et1 <- chains[[1]][,grep("^eta", colnames(chains[[1]]))]
al1 <- chains[[1]][,grep("^alpha", colnames(chains[[1]]))]
th2 <- chains[[2]][,grep("theta", colnames(chains[[2]]))]
dt2 <- chains[[2]][,grep("delta", colnames(chains[[2]]))]
bt2 <- chains[[2]][,grep("beta", colnames(chains[[2]]))]
et2 <- chains[[2]][,grep("^eta", colnames(chains[[2]]))]
al2 <- chains[[2]][,grep("^alpha", colnames(chains[[2]]))]
Theta.cjr <-  rbind (th1, th2)
Delta.cjr <-  rbind (dt1, dt2)
Beta.cjr <- rbind (bt1, bt2)
Eta.cjr <- rbind (et1, et2)
Alpha.cjr <- rbind (al1, al2)

files2load <- c("objects_2MPE_Brazil.Rdata",  
                "objects_1MPE_Brazil.Rdata")  


load (files2load[2])
Theta.P <-  items2export$theta
Delta.P <-  items2export$delta
Beta.P <- items2export$beta
Eta.P <- items2export$eta
Alpha.P <- items2export$alpha
rm (items2export)

load (files2load[1])
Theta.C <-  items2export$theta
Delta.C <-  items2export$delta
Beta.C <- items2export$beta
Eta.C <- items2export$eta
Alpha.C <- items2export$alpha

voteDisagree  <- items2export$voteDis
voteAgree     <- items2export$voteAgr
voteUncertain <- items2export$voteUnc
rm (items2export)

# PLOT ideal points
# Figure B1 (Appendix) - left
# pdf (paste0(graphPath, "/fig3.pdf"), h=7, w=7)   # formerly cjrBrazil.pdf
par (mfrow=c(1,1), mar=c(4,4,4,1))
plot ( xy.coords (colMeans(Theta.cjr)[new.party.color=="gray"]
                    , colMeans(Delta.cjr)[new.party.color=="gray"])
         , pch=1, col="gray"
       , xlim=c(-3,3), ylim=c(-3,3)
       , xlab="Left-Right", ylab="Government-Opposition"
       , main="Brazil, ideal points based on CJR model")
points (xy.coords (colMeans(Theta.cjr)[new.party.color!="gray"]
      , colMeans(Delta.cjr)[new.party.color!="gray"])
      # , col=all.party.colors
      , col=new.party.color[new.party.color!="gray"]
      #     , col=party.colors
      , pch= 19)
legend ("bottomleft", legend=c("PT","PMDB","PSDB","DEM","PP","PR","Other")
        , col=c("red","purple","black","blue","yellow","orange","grey"), lwd=5, bty="n")
# dev.off ()

# Figure B1 (Appendix) - right
# pdf (paste0(graphPath, "/fig4.pdf"), h=7, w=7)   # formerly cjrBrazilDiscrimination.pdf
par (mfrow=c(1,1), mar=c(4,4,4,1))
plot ( xy.coords (colMeans(Beta.cjr)
                  , colMeans(Eta.cjr))
       , pch=19, col="gray"
       , xlim=c(-3,3), ylim=c(-3,3)
       , xlab="Left-Right", ylab="Government-Opposition"
       , main="Discrimination parameters")
points (xy.coords (colMeans(Beta.cjr)[which (unique(Vote.Features$vote.number) %in% c("2011026","2011075","2013077","2013078"))]
                   , colMeans(Eta.cjr)[which (unique(Vote.Features$vote.number) %in% c("2011026","2011075","2013077","2013078"))])
        , pch=19
        , col="black", cex=1.5)
abline (h=c(-2,0,2), lty=3)
abline (v=c(-2,0,2), lty=3)
# dev.off ()




# The following graph does not appear either in the paper or in the appendix,
# but it reproduces the dispersion of ideal points in Brazil based on the 1MPE model
# pdf (paste0(graphPath, "/fig19.pdf"), h=7, w=7)    # formerly pemBrazil.pdf
par (mfrow=c(1,1), mar=c(4,4,4,1))
plot ( xy.coords (colMeans(Theta.P)[new.party.color=="gray"]
                  , colMeans(Delta.P)[new.party.color=="gray"])
       , pch=1, col="gray"
       , xlim=c(-3,3), ylim=c(-3,3)
       , ylab="Government-Opposition", xlab="Left-Right"
       , main="Brazil, 1MPE model")
points (xy.coords (colMeans(Theta.P)[new.party.color!="gray"]
                   , colMeans(Delta.P)[new.party.color!="gray"])
        # , col=all.party.colors
        , col=new.party.color[new.party.color!="gray"]
        #     , col=party.colors
        , pch= 19)
legend ("bottomleft", legend=c("PT","PMDB","PSDB","DEM","PP","PR","Other")
        , col=c("red","purple","black","blue","yellow","orange","grey"), lwd=5, bty="n")
# dev.off ()

# Figure H1 (appendix) - left
# pdf (paste0(graphPath, "/fig17.pdf"), h=7, w=7)   # formerly penvmBrazil.pdf
par (mfrow=c(1,1), mar=c(4,4,4,1))
plot ( xy.coords (colMeans(Theta.C)[new.party.color=="gray"]
                  , colMeans(Delta.C)[new.party.color=="gray"])
       , pch=1, col="gray"
       , xlim=c(-3,3), ylim=c(-3,3)
       , xlab="Left-Right", ylab="Government-Opposition"
       , main="Brazil, 2MPE model")
points (xy.coords (colMeans(Theta.C)[new.party.color!="gray"]
                   , colMeans(Delta.C)[new.party.color!="gray"])
        # , col=all.party.colors
        , col=new.party.color[new.party.color!="gray"]
        #     , col=party.colors
        , pch= 19)
legend ("bottomleft", legend=c("PT","PMDB","PSDB","DEM","PP","PR","Other")
        , col=c("red","purple","black","blue","yellow","orange","grey"), lwd=5, bty="n")
# dev.off ()

## Alternative ideal point plots
dat1 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="red"]
                   , y = colMeans(Delta.cjr)[new.party.color=="red"])
dat2 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="black"]
                   , y = colMeans(Delta.cjr)[new.party.color=="black"])
dat3 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="orange"]
                   , y = colMeans(Delta.cjr)[new.party.color=="orange"])
dat4 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="purple"]
                   , y = colMeans(Delta.cjr)[new.party.color=="purple"])
dat5 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="blue"]
                   , y = colMeans(Delta.cjr)[new.party.color=="blue"])
dat6 <- data.frame(x = colMeans(Theta.cjr)[new.party.color=="yellow"]
                   , y = colMeans(Delta.cjr)[new.party.color=="yellow"])
dat <- cbind(rbind(dat1, dat2, dat3, dat4, dat5, dat6)
             , df = rep.int(c("PT","PSDB","PR","PMDB","DEM","PP")
                  , c(nrow(dat1), nrow(dat2), nrow(dat3)
                      , nrow(dat4), nrow(dat5), nrow(dat6))))




datP1 <- data.frame(x = colMeans(Theta.P)[new.party.color=="red"]
                   , y = colMeans(Delta.P)[new.party.color=="red"])
datP2 <- data.frame(x = colMeans(Theta.P)[new.party.color=="black"]
                   , y = colMeans(Delta.P)[new.party.color=="black"])
datP3 <- data.frame(x = colMeans(Theta.P)[new.party.color=="orange"]
                   , y = colMeans(Delta.P)[new.party.color=="orange"])
datP4 <- data.frame(x = colMeans(Theta.P)[new.party.color=="purple"]
                   , y = colMeans(Delta.P)[new.party.color=="purple"])
datP5 <- data.frame(x = colMeans(Theta.P)[new.party.color=="blue"]
                   , y = colMeans(Delta.P)[new.party.color=="blue"])
datP6 <- data.frame(x = colMeans(Theta.P)[new.party.color=="yellow"]
                   , y = colMeans(Delta.P)[new.party.color=="yellow"])

datC1 <- data.frame(x = colMeans(Theta.C)[new.party.color=="red"]
                    , y = colMeans(Delta.C)[new.party.color=="red"])
datC2 <- data.frame(x = colMeans(Theta.C)[new.party.color=="black"]
                    , y = colMeans(Delta.C)[new.party.color=="black"])
datC3 <- data.frame(x = colMeans(Theta.C)[new.party.color=="orange"]
                    , y = colMeans(Delta.C)[new.party.color=="orange"])
datC4 <- data.frame(x = colMeans(Theta.C)[new.party.color=="purple"]
                    , y = colMeans(Delta.C)[new.party.color=="purple"])
datC5 <- data.frame(x = colMeans(Theta.C)[new.party.color=="blue"]
                    , y = colMeans(Delta.C)[new.party.color=="blue"])
datC6 <- data.frame(x = colMeans(Theta.C)[new.party.color=="yellow"]
                    , y = colMeans(Delta.C)[new.party.color=="yellow"])


require(mixtools)
med.th.1 <- med.dl.1 <- c()
med.th.2 <- med.dl.2 <- c()
med.th.3 <- med.dl.3 <- c()
med.th.4 <- med.dl.4 <- c()
med.th.5 <- med.dl.5 <- c()
med.th.6 <- med.dl.6 <- c()
for (i in 1:nrow(Theta.cjr))  {
  med.th.1[i] <- median(Theta.cjr[i,new.party.color=="red"])
  med.dl.1[i] <- median(Delta.cjr[i,new.party.color=="red"])
  med.th.2[i] <- median(Theta.cjr[i,new.party.color=="black"])
  med.dl.2[i] <- median(Delta.cjr[i,new.party.color=="black"])
  med.th.3[i] <- median(Theta.cjr[i,new.party.color=="orange"])
  med.dl.3[i] <- median(Delta.cjr[i,new.party.color=="orange"])
  med.th.4[i] <- median(Theta.cjr[i,new.party.color=="purple"])
  med.dl.4[i] <- median(Delta.cjr[i,new.party.color=="purple"])
  med.th.5[i] <- median(Theta.cjr[i,new.party.color=="blue"])
  med.dl.5[i] <- median(Delta.cjr[i,new.party.color=="blue"])
  med.th.6[i] <- median(Theta.cjr[i,new.party.color=="yellow"])
  med.dl.6[i] <- median(Delta.cjr[i,new.party.color=="yellow"])
}

med.th.C1 <- med.dl.C1 <- c()
med.th.C2 <- med.dl.C2 <- c()
med.th.C3 <- med.dl.C3 <- c()
med.th.C4 <- med.dl.C4 <- c()
med.th.C5 <- med.dl.C5 <- c()
med.th.C6 <- med.dl.C6 <- c()
for (i in 1:nrow(Theta.C))  {
  med.th.C1[i] <- median(Theta.C[i,new.party.color=="red"])
  med.dl.C1[i] <- median(Delta.C[i,new.party.color=="red"])
  med.th.C2[i] <- median(Theta.C[i,new.party.color=="black"])
  med.dl.C2[i] <- median(Delta.C[i,new.party.color=="black"])
  med.th.C3[i] <- median(Theta.C[i,new.party.color=="orange"])
  med.dl.C3[i] <- median(Delta.C[i,new.party.color=="orange"])
  med.th.C4[i] <- median(Theta.C[i,new.party.color=="purple"])
  med.dl.C4[i] <- median(Delta.C[i,new.party.color=="purple"])
  med.th.C5[i] <- median(Theta.C[i,new.party.color=="blue"])
  med.dl.C5[i] <- median(Delta.C[i,new.party.color=="blue"])
  med.th.C6[i] <- median(Theta.C[i,new.party.color=="yellow"])
  med.dl.C6[i] <- median(Delta.C[i,new.party.color=="yellow"])
}


dat1.med <- data.frame(x = med.th.1, y = med.dl.1)
dat2.med <- data.frame(x = med.th.2, y = med.dl.2)
dat3.med <- data.frame(x = med.th.3, y = med.dl.3)
dat4.med <- data.frame(x = med.th.4, y = med.dl.4)
dat5.med <- data.frame(x = med.th.5, y = med.dl.5)
dat6.med <- data.frame(x = med.th.6, y = med.dl.6)

dat1.med.C <- data.frame(x = med.th.C1, y = med.dl.C1)
dat2.med.C <- data.frame(x = med.th.C2, y = med.dl.C2)
dat3.med.C <- data.frame(x = med.th.C3, y = med.dl.C3)
dat4.med.C <- data.frame(x = med.th.C4, y = med.dl.C4)
dat5.med.C <- data.frame(x = med.th.C5, y = med.dl.C5)
dat6.med.C <- data.frame(x = med.th.C6, y = med.dl.C6)

# Find party medians, to draw convex hull
Brazil.party.medians <- cbind (c(median (dat1$x), median (dat2$x), median (dat3$x), median (dat4$x), median (dat5$x), median (dat6$x)),
       c(median (dat1$y), median (dat2$y), median (dat3$y), median (dat4$y), median (dat5$y), median (dat6$y)))

Brazil.party.medians.C <- cbind (c(median (datC1$x), median (datC2$x), median (datC3$x), median (datC4$x), median (datC5$x), median (datC6$x)),
       c(median (datC1$y), median (datC2$y), median (datC3$y), median (datC4$y), median (datC5$y), median (datC6$y)))


# Calculate within-party distances to mean position

sd.CJR <- c(sd (sqrt ((dat1$x - mean (dat1$x))^2 + (dat1$y - mean (dat1$y))^2))
                  , sd (sqrt ((dat2$x - mean (dat2$x))^2 + (dat2$y - mean (dat2$y))^2))
                  , sd (sqrt ((dat3$x - mean (dat3$x))^2 + (dat3$y - mean (dat3$y))^2))
                  , sd (sqrt ((dat4$x - mean (dat4$x))^2 + (dat4$y - mean (dat4$y))^2))
                  , sd (sqrt ((dat5$x - mean (dat5$x))^2 + (dat5$y - mean (dat5$y))^2))
                  , sd (sqrt ((dat6$x - mean (dat6$x))^2 + (dat6$y - mean (dat6$y))^2)))

sd.PENUM <- c(sd (sqrt ((datC1$x - mean (datC1$x))^2 + (datC1$y - mean (datC1$y))^2))
                  , sd (sqrt ((datC2$x - mean (datC2$x))^2 + (datC2$y - mean (datC2$y))^2))
                  , sd (sqrt ((datC3$x - mean (datC3$x))^2 + (datC3$y - mean (datC3$y))^2))
                  , sd (sqrt ((datC4$x - mean (datC4$x))^2 + (datC4$y - mean (datC4$y))^2))
                  , sd (sqrt ((datC5$x - mean (datC5$x))^2 + (datC5$y - mean (datC5$y))^2))
                  , sd (sqrt ((datC6$x - mean (datC6$x))^2 + (datC6$y - mean (datC6$y))^2)))

sd.Brazil <- as.table (rbind (sd.CJR, sd.PENUM))
colnames (sd.Brazil) <- c("PT","PSDB","PR","PMDB","DEM","PP")
print (sd.Brazil)


# Function to plot convex hull, from https://chitchatr.wordpress.com/2011/12/30/convex-hull-around-scatter-plot-in-r/
Plot_ConvexHull<-function(xcoord, ycoord, lcolor, lwidth){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor, lwd = lwidth)
}  

Area_ConvexHull <- function (coordMatrix) {
  require (sp)
  box.hpts <- chull(x = coordMatrix[,1], y = coordMatrix[,2])
  box.hpts <- c(box.hpts, box.hpts[1])
  box.chull.coords <- coordMatrix[box.hpts,]
  chull.poly <- Polygon(box.chull.coords, hole=F)
  chull.area <- chull.poly@area 
  return (chull.area)
}


# Fig 1
# pdf (paste0 (graphPath, "/fig1.pdf"), h=7, w=7)  #formerly BrazilPartyMediansConvexHull.pdf
par (las=0, mar=c(3,3,0,0))
plot (colMeans(Theta.cjr), colMeans(Delta.cjr)
      , type="n"
      #     , main="2 bill spike 2 bill semispike"
      , xlim=c(-2.25,1.25), ylim=c(-1.5,2)
      , ylab="", xlab="", axes=F)
axis (1, at=c(-2.5,2.5)
      , labels=c(NA,NA))
axis (2, at=c(-2.5,2.5)
      , labels=c(NA,NA))
mtext (side=1, text="Left-Right", line=1)
mtext (side=2, text="Government-Opposition", line=1)
Plot_ConvexHull(xcoord = Brazil.party.medians[,1], ycoord = Brazil.party.medians[,2], lcolor = "gray", lwidth = 3)
Plot_ConvexHull(xcoord = Brazil.party.medians.C[,1], ycoord = Brazil.party.medians.C[,2], lcolor = "black", lwidth = 3)
points (xy.coords(Brazil.party.medians), pch=19, col="gray", cex=1.2)
points (xy.coords(Brazil.party.medians.C), pch=19, col="black", cex=1.2)
text (xy.coords(Brazil.party.medians), labels=c("PT","PSDB","PR","PMDB","DEM","PP")
      , pos=c(3,1,3,3,4,2), col="dimgray")
text (xy.coords(Brazil.party.medians.C), labels=c("PT","PSDB","PR","PMDB","DEM","PP")
      , pos=c(2,2,2,4,2,4), col="black")
legend ("bottomleft"
        , legend = c(paste0("Party effects convex hull area: ", round (Area_ConvexHull (Brazil.party.medians.C), 2))
                  , paste0("IRT convex hull area: ", round (Area_ConvexHull (Brazil.party.medians), 2)))
        , bty="n", lwd=3, col=c("black","gray"))
# dev.off()


## Table of convex hull areas for parties
ConvexHull.Brazil <- rbind (c(Area_ConvexHull (dat1), Area_ConvexHull (dat2), Area_ConvexHull (dat3)
                              , Area_ConvexHull (dat4), Area_ConvexHull (dat5), Area_ConvexHull (dat6)
                              , mean (c(Area_ConvexHull (dat1), Area_ConvexHull (dat2), Area_ConvexHull (dat3)
                                        , Area_ConvexHull (dat4), Area_ConvexHull (dat5), Area_ConvexHull (dat6))))
                            , c(Area_ConvexHull (datC1), Area_ConvexHull (datC2), Area_ConvexHull (datC3)
                                , Area_ConvexHull (datC4), Area_ConvexHull (datC5), Area_ConvexHull (datC6)
                                , mean(c(Area_ConvexHull (datC1), Area_ConvexHull (datC2), Area_ConvexHull (datC3)
                                         , Area_ConvexHull (datC4), Area_ConvexHull (datC5), Area_ConvexHull (datC6)))))
colnames (ConvexHull.Brazil) <- c("PT","PSDB","PR","PMDB","DEM","PP","Avg.")
rownames (ConvexHull.Brazil) <- c("IRT","Party Effects")


# Table 7, Brazil half
xtable (ConvexHull.Brazil, align=c("l","c","c","c","c","c","c","c")
        , caption="Convex hull area of Brazilian parties based on IRT and Party Effects models"
        , label="T:convexHullBrazil")






###################
### Plot Mexico ###
###################
rm (list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)


# Set graph path
graphPath <- getwd()

##plot CJR model
load(file="mexico_CJR_auxiliary_output.Rda")
chains <- chainsCJR  
th1 <- chains[[1]][,grep("theta", colnames(chains[[1]]))]
dt1 <- chains[[1]][,grep("delta", colnames(chains[[1]]))]
bt1 <- chains[[1]][,grep("beta", colnames(chains[[1]]))]
et1 <- chains[[1]][,grep("^eta", colnames(chains[[1]]))]
al1 <- chains[[1]][,grep("^alpha", colnames(chains[[1]]))]
th2 <- chains[[2]][,grep("theta", colnames(chains[[2]]))]
dt2 <- chains[[2]][,grep("delta", colnames(chains[[2]]))]
bt2 <- chains[[2]][,grep("beta", colnames(chains[[2]]))]
et2 <- chains[[2]][,grep("^eta", colnames(chains[[2]]))]
al2 <- chains[[2]][,grep("^alpha", colnames(chains[[2]]))]
Theta.cjr <-  rbind (th1, th2)
Delta.cjr <-  rbind (dt1, dt2)
Beta.cjr <- rbind (bt1, bt2)
Eta.cjr <- rbind (et1, et2)
Alpha.cjr <- rbind (al1, al2)


##plot Table F7
load ("objects_1MPE_Mexico.Rdata")
Theta.P <-  items2export$theta
Delta.P <-  items2export$delta
Beta.P  <-  items2export$beta
Eta.P   <-  items2export$eta
Alpha.P <-  items2export$alpha
rm (items2export)




load ("objects_2MPE_Mexico.Rdata")
Theta.C <-  items2export$theta
Delta.C <-  items2export$delta
Beta.C  <-  items2export$beta
Eta.C   <-  items2export$eta
Alpha.C <-  items2export$alpha
rm (items2export)




load ("mexico_rollcall_long_format_mean.RData")

# PLOT ideal points
party.color <- as.numeric (by (good.comp.party$party, good.comp.party$final.leg.no, unique))
party.color <- car::recode (party.color, "1='green';
                            2='blue';
                            3='red';
                            4='yellow';
                            5='orange';
                            6='black';
                            7='purple';
                            8='gray'")



dat1 <- data.frame(x = colMeans(Theta.cjr)[party.color=="red"]
                   , y = colMeans(Delta.cjr)[party.color=="red"])
dat2 <- data.frame(x = colMeans(Theta.cjr)[party.color=="blue"]
                   , y = colMeans(Delta.cjr)[party.color=="blue"])
dat3 <- data.frame(x = colMeans(Theta.cjr)[party.color=="green"]
                   , y = colMeans(Delta.cjr)[party.color=="green"])

datP1 <- data.frame(x = colMeans(Theta.P)[party.color=="red"]
                    , y = colMeans(Delta.P)[party.color=="red"])
datP2 <- data.frame(x = colMeans(Theta.P)[party.color=="blue"]
                    , y = colMeans(Delta.P)[party.color=="blue"])
datP3 <- data.frame(x = colMeans(Theta.P)[party.color=="green"]
                    , y = colMeans(Delta.P)[party.color=="green"])

datC1 <- data.frame(x = colMeans(Theta.C)[party.color=="red"]
                    , y = colMeans(Delta.C)[party.color=="red"])
datC2 <- data.frame(x = colMeans(Theta.C)[party.color=="blue"]
                    , y = colMeans(Delta.C)[party.color=="blue"])
datC3 <- data.frame(x = colMeans(Theta.C)[party.color=="green"]
                    , y = colMeans(Delta.C)[party.color=="green"])


party.names <- c("PRI","PAN","PRD","PVEM","PANAL","MC","PT","Indep.")

# Figure B2 (Appendix) left
# pdf (paste0 (graphPath, "/fig5.pdf"), h=7, w=7)  # formerly cjrMexico.pdf
par (mar=c(4,4,4,1), las=1)
plot (colMeans(Theta.cjr), colMeans(Delta.cjr)
      , col=party.color
      , pch=19
      , main="Mexico, ideal points based on CJR model"
      , xlim=c(-4,4), ylim=c(-4,4)
      , ylab="", xlab="")
par (las=0)
mtext (1, line=2.5, text="Economic-distributive divide")
mtext (2, line=2.5, text="Political organization")
legend ("topleft", legend=party.names
        , col=c("green","blue","red",
                "yellow","orange","black","purple",
                "gray"), pch=19, bty="n", cex=0.8)
# dev.off()

# The following graph does not appear either in the paper or in the appendix,
# but it reproduces the dispersion of ideal points in Mexico based on the 1MPE model
# pdf (paste0 (graphPath, "/fig20.pdf"), h=7, w=7)  
par (mfrow=c(1,1))
plot (colMeans(Theta.P), colMeans(Delta.P)
      , col=party.color
      , pch=19
      , main="Mexico, 1MPE model"
      , xlim=c(-4,4), ylim=c(-4,4)
      , ylab="", xlab="")
par (las=0)
mtext (1, line=2.5, text="Economic-distributive divide")
mtext (2, line=2.5, text="Political organization")
legend ("topleft", legend=party.names
        , col=c("green","blue","red",
                "yellow","orange","black","purple",
                "gray"), pch=19, bty="n", cex=0.8)
# dev.off ()


# Figure H1 (Appendix) Right/Mexico
# pdf (paste0 (graphPath, "/fig18.pdf"), h=7, w=7)  
par (mfrow=c(1,1))
plot (colMeans(Theta.C), colMeans(Delta.C)
      , col=party.color
      , pch=19
      , main="Mexico, 2MPE model"
      , xlim=c(-4,4), ylim=c(-4,4)
      , ylab="", xlab="")
par (las=0)
mtext (1, line=2.5, text="Economic-distributive divide")
mtext (2, line=2.5, text="Political organization")
legend ("topleft", legend=party.names
        , col=c("green","blue","red",
        "yellow","orange","black","purple",
        "gray"), pch=19, bty="n", cex=0.8)
# dev.off ()


# Figure B2 (Appendix) Right
# pdf (paste0 (graphPath, "/fig6.pdf"), h=7, w=7)  # formerly cjrMexicoItemParams.pdf
par (mfrow=c(1,1), mar=c(4,4,4,1))
plot (colMeans(Beta.cjr), colMeans(Eta.cjr)
      #     , col=all.party.colors
      # , col=party.color
      #, col=party.colors
      , pch=19, col="gray"
      , main="Discrimination parameters"
      , xlim=c(-4,4), ylim=c(-4,4)
      , ylab="", xlab="")
# points (xy.coords (colMeans(Beta.cjr)[which (apply (Beta.cjr, 2, sd)==0)]
#         , colMeans(Eta.cjr)[which (apply (Beta.cjr, 2, sd)==0)]), pch=19)
par (las=0)
mtext (1, line=2.5, text="Economic-distributive divide")
mtext (2, line=2.5, text="Political organization")
points (xy.coords (colMeans(Beta.cjr)[c(89,150,174)]     
                   , colMeans(Eta.cjr)[c(89,150,174)])       
        , pch=19
        , col="black", cex=1.5)
abline (h=c(-2,0,2), lty=3)
abline (v=c(-2,0,2), lty=3)
# dev.off()




# Mexico medians
med.th.1 <- med.dl.1 <- c()
med.th.2 <- med.dl.2 <- c()
med.th.3 <- med.dl.3 <- c()
for (i in 1:nrow(Theta.cjr))  {
  med.th.1[i] <- median(Theta.cjr[i,party.color=="red"])
  med.dl.1[i] <- median(Delta.cjr[i,party.color=="red"])
  med.th.2[i] <- median(Theta.cjr[i,party.color=="blue"])
  med.dl.2[i] <- median(Delta.cjr[i,party.color=="blue"])
  med.th.3[i] <- median(Theta.cjr[i,party.color=="green"])
  med.dl.3[i] <- median(Delta.cjr[i,party.color=="green"])
}

med.th.C1 <- med.dl.C1 <- c()
med.th.C2 <- med.dl.C2 <- c()
med.th.C3 <- med.dl.C3 <- c()
for (i in 1:nrow(Theta.C))  {
  med.th.C1[i] <- median(Theta.C[i,party.color=="red"])
  med.dl.C1[i] <- median(Delta.C[i,party.color=="red"])
  med.th.C2[i] <- median(Theta.C[i,party.color=="blue"])
  med.dl.C2[i] <- median(Delta.C[i,party.color=="blue"])
  med.th.C3[i] <- median(Theta.C[i,party.color=="green"])
  med.dl.C3[i] <- median(Delta.C[i,party.color=="green"])
}


dat1.med <- data.frame(x = med.th.1, y = med.dl.1)
dat2.med <- data.frame(x = med.th.2, y = med.dl.2)
dat3.med <- data.frame(x = med.th.3, y = med.dl.3)

dat1.med.C <- data.frame(x = med.th.C1, y = med.dl.C1)
dat2.med.C <- data.frame(x = med.th.C2, y = med.dl.C2)
dat3.med.C <- data.frame(x = med.th.C3, y = med.dl.C3)


# Find party medians, to draw convex hull
Mexico.party.medians <- cbind (c(median (dat1$x), median (dat2$x), median (dat3$x)),
                               c(median (dat1$y), median (dat2$y), median (dat3$y)))

Mexico.party.medians.C <- cbind (c(median (datC1$x), median (datC2$x), median (datC3$x)),
                                 c(median (datC1$y), median (datC2$y), median (datC3$y)))

# Calculate within-party distances to mean position

sd.CJR <- c(sd (sqrt ((dat1$x - mean (dat1$x))^2 + (dat1$y - mean (dat1$y))^2))
            , sd (sqrt ((dat2$x - mean (dat2$x))^2 + (dat2$y - mean (dat2$y))^2))
            , sd (sqrt ((dat3$x - mean (dat3$x))^2 + (dat3$y - mean (dat3$y))^2)))

sd.PENUM <- c(sd (sqrt ((datC1$x - mean (datC1$x))^2 + (datC1$y - mean (datC1$y))^2))
              , sd (sqrt ((datC2$x - mean (datC2$x))^2 + (datC2$y - mean (datC2$y))^2))
              , sd (sqrt ((datC3$x - mean (datC3$x))^2 + (datC3$y - mean (datC3$y))^2)))

sd.Mexico <- as.table (rbind (sd.CJR, sd.PENUM))
colnames (sd.Mexico) <- c("PRD","PAN","PRI")



# Function to plot convex hull, from https://chitchatr.wordpress.com/2011/12/30/convex-hull-around-scatter-plot-in-r/
Plot_ConvexHull<-function(xcoord, ycoord, lcolor, lwidth){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor, lwd = lwidth)
}  

Area_ConvexHull <- function (coordMatrix) {
  require (sp)
  box.hpts <- chull(x = coordMatrix[,1], y = coordMatrix[,2])
  box.hpts <- c(box.hpts, box.hpts[1])
  box.chull.coords <- coordMatrix[box.hpts,]
  chull.poly <- Polygon(box.chull.coords, hole=F)
  chull.area <- chull.poly@area 
  return (chull.area)
}


# Figure 1, main paper, right
# pdf (paste0(graphPath, "/fig2.pdf"), h=7, w=7)  # formerly MexicoPartyMediansConvexHull.pdf
par (las=0, mar=c(3,3,0,0))
plot (colMeans(Theta.cjr), colMeans(Delta.cjr)
      , type="n"
      #      , main="2 bill spike 2 bill semispike"
      , xlim=c(-1.75,1.75), ylim=c(-2,1.5)
      , ylab="", xlab="", axes=F)
axis (1, at=c(-2.5,2.5)
      , labels=c(NA,NA))
axis (2, at=c(-2.5,2.5)
      , labels=c(NA,NA))
mtext (side=1, text="Economic-Distributive Divide", line=1)
mtext (side=2, text="Political Organization", line=1)
Plot_ConvexHull(xcoord = Mexico.party.medians[,1], ycoord = Mexico.party.medians[,2], lcolor = "gray", lwidth = 3)
Plot_ConvexHull(xcoord = Mexico.party.medians.C[,1], ycoord = Mexico.party.medians.C[,2], lcolor = "black", lwidth = 3)
points (xy.coords(Mexico.party.medians), pch=19, col="gray", cex=1.2)
points (xy.coords(Mexico.party.medians.C), pch=19, col="black", cex=1.2)
text (xy.coords(Mexico.party.medians), labels=c("PRD","PAN","PRI")
      , pos=c(2,1,3), col="gray")
text (xy.coords(Mexico.party.medians.C), labels=c("PRD","PAN","PRI")
      , pos=c(2,1,3), col="black")
legend ("bottomleft"
        , legend = c(paste0("Party effects convex hull area: ", round (Area_ConvexHull (Mexico.party.medians.C), 1))
                     , paste0("IRT convex hull area: ", round (Area_ConvexHull (Mexico.party.medians), 2)))
        , bty="n", lwd=3, col=c("black","gray"))
# dev.off()




## Table of convex hull areas for parties
ConvexHull.Mexico <- rbind (c(Area_ConvexHull (dat1), Area_ConvexHull (dat2), Area_ConvexHull (dat3),
                              mean (c(Area_ConvexHull (dat1), Area_ConvexHull (dat2), Area_ConvexHull (dat3)))),
                            c(Area_ConvexHull (datC1), Area_ConvexHull (datC2), Area_ConvexHull (datC3),
                              mean(c(Area_ConvexHull (datC1), Area_ConvexHull (datC2), Area_ConvexHull (datC3)))))
colnames (ConvexHull.Mexico) <- c("PRD","PAN","PRI","Avg.")
rownames (ConvexHull.Mexico) <- c("IRT","Party Effects")


# Table 7, Mexico half
xtable (ConvexHull.Mexico, align=c("l","c","c","c","c")
        , caption="Convex hull area of Mexican parties based on IRT and Party Effects models"
        , label="T:convexHullMexico")   # should be T:cohesive





##Descriptive statistics
rm (list=ls())

load ("brazil_rollcall_long_format_mean.Rdata")

Procedure.leader.pressure<- Procedure.PT+Procedure.PSDB+Procedure.DEM+Procedure.PMDB+Procedure.PP+Procedure.PR
Procedure.leader<- pmax(Vote.Features$PT.leader.request.procedure.x, Vote.Features$PSDB.leader.request.procedure.x, Vote.Features$DEM.leader.request.procedure.x, Vote.Features$PMDB.leader.request.procedure.x, Vote.Features$PP.leader.request.procedure.x, Vote.Features$PR.leader.request.procedure.x)

pcb <- ifelse(Vote.Features$pcb.request.x==1, 1, 0) 


P8=seniority
P10=committee.chair
P12=district.income
P14=percentile.rank
P16=log.district.M
P18=molten.roll.call$year
P20=Procedure.leader
P22=minister
P24=Vote.Features$president.position.x
P26=pcb


var.avg<- rbind(P8, P10, P12, P14, P16, P18, P20, P22, P24, P26)
var<- round(t(apply(var.avg, 1, summary)),3)
rownames(var)<- c("Seniority", "Committee Chair", "District Income", "Percentile Rank", "District Magnitude", "Year", "Procedural Vote", "Cabinet Posts", "President Position", "Party Agreements")


#  Table D1 (Appendix) - Brazil
xtable(var[,c(1,4,6)] 
       , caption="Descriptive statistics for variables included in Table D1 (Brazil)"
       , label="T:Descriptive")



#  Mexico
rm (list=ls())

load ("mexico_rollcall_long_format_mean.Rdata")

P5=good.R67$pressure   #PR
P7=good.R64$pressure  #poverty dummy
P9=good.R65$pressure    # year
P11=good.R75$pressure   #victory margin dummy
P13=good.R72$pressure  #sta.muni.leg
P15=good.R77$pressure  #CabinetExp
P17=good.R71$pressure   #prezposition
P19=good.R73$pressure    #commPres
P21=good.R56$pressure    #pcbrequest


var.avg<- rbind(P5, P7, P9, P11, P13, P15, P17, P19, P21)
var<- round(t(apply(var.avg, 1, summary)),3)
rownames(var)<- c("PR", "Poverty", "Year", "Victory Margin", "State Legislators", "Cabinet Posts", "President Position", "Committee Chair" , "Party Agreements")

#  Table D1 (Appendix) Mexico
xtable(var[,c(1,4,6)] 
       , caption="Descriptive statistics for variables included in Tables D1 (Mexico)"
       , label="T:Descriptive")



##################################################
#### Code to produce Table C2 in the appendix ####
##################################################


#### MEXICO ####

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()

# Vectorized information on vote, party, and bill
load ("mexico_rollcall_long_format_mean.Rdata") 
vectorVotes <- data.frame (cbind (good.comp.party[,1:3], good.RC$RC))
vectorVotes$vote <- as.numeric (vectorVotes$vote)

# vectorVotes contains a vectorized version of legislator-vote information
rm (list=setdiff(ls(),
                 c("vectorVotes")))

# Use vectorVotes to reconstruct party unity indexes
unity.index.function <- function (tmp) {
  no.Aye <- length (tmp[tmp==1 & !is.na(tmp)])
  no.Nay <- length (tmp[tmp==0 & !is.na(tmp)])
  no.Abs <- length (tmp[is.na(tmp)])
  rice.wo.abs <- abs (no.Aye-no.Nay) / (no.Aye+no.Nay) # Rice excluding abstentions (ties are 0)
  rice.w.abs  <- abs (no.Aye-no.Nay) / (no.Aye+no.Nay+no.Abs) # Rice including abstentions
  alt.rice.wo.abs <- max (no.Aye, no.Nay) / (no.Aye+no.Nay) # Ties are 0.5
  alt.rice.w.abs  <- max (no.Aye, no.Nay) / (no.Aye+no.Nay+no.Abs)
  agreement.index <- (max (no.Aye, no.Nay, no.Abs) - 0.5*(length(tmp) - max (no.Aye, no.Nay, no.Abs))) / length(tmp)   
  absent <- no.Abs/(no.Aye + no.Nay + no.Abs)
  all.unity.indexes <- c (rice.wo.abs, rice.w.abs, alt.rice.wo.abs, alt.rice.w.abs, agreement.index, absent)
  return (all.unity.indexes)
}


PRI.unity.index.vector <- PAN.unity.index.vector <- PRD.unity.index.vector <- c()
for (i in 1:length(unique(vectorVotes$vote))){
  tmp.vote.pri <- vectorVotes$good.RC.RC[vectorVotes$vote==i & vectorVotes$party==1]
  tmp.vote.pan <- vectorVotes$good.RC.RC[vectorVotes$vote==i & vectorVotes$party==2]
  tmp.vote.prd <- vectorVotes$good.RC.RC[vectorVotes$vote==i & vectorVotes$party==3]
  tmp.col.pri  <- unity.index.function (tmp.vote.pri)
  tmp.col.pan  <- unity.index.function (tmp.vote.pan)
  tmp.col.prd  <- unity.index.function (tmp.vote.prd)
  PRI.unity.index.vector <- cbind (PRI.unity.index.vector, tmp.col.pri)
  PAN.unity.index.vector <- cbind (PAN.unity.index.vector, tmp.col.pan)
  PRD.unity.index.vector <- cbind (PRD.unity.index.vector, tmp.col.prd)
}

rownames (PRI.unity.index.vector) <- rownames (PAN.unity.index.vector) <- rownames (PRD.unity.index.vector) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index","absent.pct")

round(mean(PRI.unity.index.vector[6,]), 2); round(median(PRI.unity.index.vector[6,]), 2)
round(mean(PAN.unity.index.vector[6,]), 2); round(median(PAN.unity.index.vector[6,]), 2)
round(mean(PRD.unity.index.vector[6,]), 2); round(median(PRD.unity.index.vector[6,]), 2)


# Load Party Pressure object for Mexico
load("mexico_ptirt_mcmc.Rda")
X.e   <- rbind (chainsPTIRT[[1]][,grep("X",  colnames (chainsPTIRT[[1]]))],
                chainsPTIRT[[2]][,grep("X",  colnames (chainsPTIRT[[2]]))])

X.pri  <- colMeans (X.e[,grep(",1]", colnames(X.e))])
X.pan  <- colMeans (X.e[,grep(",2]", colnames(X.e))])
X.prd  <- colMeans (X.e[,grep(",3]", colnames(X.e))])
rm (X.e)

party.pressure.correlations.MX <- cbind (round (cor (t(rbind (X.pri, PRI.unity.index.vector)), method="spearman")[,1], 2)
                                       , round (cor (t(rbind (X.pan, PAN.unity.index.vector)), method="spearman")[,1], 2)
                                       , round (cor (t(rbind (X.prd, PRD.unity.index.vector)), method="spearman")[,1], 2))
colnames (party.pressure.correlations.MX) <- c("PRI","PAN","PRD")
rownames (party.pressure.correlations.MX) <- c("party.pressure",
                                               "rice.wo.abs",
                                               "rice.w.abs",
                                               "alt.rice.wo.abs",
                                               "alt.rice.w.abs",
                                               "agreement.index",
                                               "absent.pct")


#### BRAZIL ####

# Load Party Pressure object for  Brazil
load ("brazil_ptirt_mcmc.Rda") 

X.e   <- rbind (chainsPTIRT[[1]][,grep("X",  colnames (chainsPTIRT[[1]]))],
                chainsPTIRT[[2]][,grep("X",  colnames (chainsPTIRT[[2]]))])

X.pt    <- X.e[,grep(",1]", colnames(X.e))]
X.pmdb  <- X.e[,grep(",2]", colnames(X.e))]
X.psdb  <- X.e[,grep(",3]", colnames(X.e))]
X.dem   <- X.e[,grep(",4]", colnames(X.e))]
X.pp    <- X.e[,grep(",5]", colnames(X.e))]
X.pr    <- X.e[,grep(",6]", colnames(X.e))]
rm (X.e)

PT.pressure <- colMeans (X.pt)
PMDB.pressure <- colMeans (X.pmdb)
PSDB.pressure <- colMeans (X.psdb)
DEM.pressure <- colMeans (X.dem)
PP.pressure <- colMeans (X.pp)
PR.pressure <- colMeans (X.pr)

# Load Brazilian party unity indexes
load ("prepared_roll_call_data_brazil.Rdata")

rownames (PT.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")
rownames (PMDB.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")
rownames (PSDB.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")
rownames (DEM.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")
rownames (PP.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")
rownames (PR.unity.index) <- c("rice.wo.abs","rice.w.abs","alt.rice.wo.abs","alt.rice.w.abs","agreement.index")

PT.scores <- rbind (PT.unity.index, PT.pressure)
PMDB.scores <- rbind (PMDB.unity.index, PMDB.pressure)
PSDB.scores <- rbind (PSDB.unity.index, PSDB.pressure)
DEM.scores <- rbind (DEM.unity.index, DEM.pressure)
PP.scores <- rbind (PP.unity.index, PP.pressure)
PR.scores <- rbind (PR.unity.index, PR.pressure)
rownames (PMDB.scores)[6] <- rownames (PSDB.scores)[6] <- rownames (DEM.scores)[6] <- rownames (PP.scores)[6] <- rownames (PR.scores)[6] <- rownames (PT.scores)[6] <- "pressure"

party.pressure.correlations.BZ <- cbind (
  round (cor (t(PT.scores), use="p", method="spearman"), 2)[-6,6],
  round (cor (t(PMDB.scores), use="p", method="spearman"), 2)[-6,6],
  round (cor (t(PSDB.scores), use="p", method="spearman"), 2)[-6,6],
  round (cor (t(DEM.scores), use="p", method="spearman"), 2)[-6,6],
  round (cor (t(PP.scores), use="p", method="spearman"), 2)[-6,6],
  round (cor (t(PR.scores), use="p", method="spearman"), 2)[-6,6]
)
colnames (party.pressure.correlations.BZ) <- c("PT","PMDB","PSDB","DEM","PP","PR")


## Table C2 (Appendix)
xtable (cbind (party.pressure.correlations.BZ, party.pressure.correlations.MX[-c(1,7),])
        , caption="Correlation between party-bill-specific pressure and unity scores for six Brazilian and three Mexican political parties (Spearman rank-based correlation indices)"
        , label="t:correlationRice"
        , digits=2)




#############################################
#### Code for Tables 6 and 8, main paper ####
#############################################

library (xtable)
library (sp)
library (coda)
library (Ternary)

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()

##plot CJR model
load(file="mexico_CJR_auxiliary_output.Rda")
chains <- chainsCJR  

th1 <- chains[[1]][,grep("theta", colnames(chains[[1]]))]
dt1 <- chains[[1]][,grep("delta", colnames(chains[[1]]))]
bt1 <- chains[[1]][,grep("beta", colnames(chains[[1]]))]
et1 <- chains[[1]][,grep("^eta", colnames(chains[[1]]))]
al1 <- chains[[1]][,grep("^alpha", colnames(chains[[1]]))]
dev1 <- chains[[1]][,grep("deviance", colnames(chains[[1]]))]
th2 <- chains[[2]][,grep("theta", colnames(chains[[2]]))]
dt2 <- chains[[2]][,grep("delta", colnames(chains[[2]]))]
bt2 <- chains[[2]][,grep("beta", colnames(chains[[2]]))]
et2 <- chains[[2]][,grep("^eta", colnames(chains[[2]]))]
al2 <- chains[[2]][,grep("^alpha", colnames(chains[[2]]))]
dev2 <- chains[[2]][,grep("deviance", colnames(chains[[2]]))]
Theta.cjr <-  rbind (th1, th2)
Delta.cjr <-  rbind (dt1, dt2)
Beta.cjr <- rbind (bt1, bt2)
Eta.cjr <- rbind (et1, et2)
Alpha.cjr <- rbind (al1, al2)
Deviance.cjr <- c (dev1, dev2)
rm (chainsCJR)

load ("objects_1MPE_Mexico.Rdata")
Theta.P <-  items2export$theta
Delta.P <-  items2export$delta
Beta.P  <-  items2export$beta
Eta.P   <-  items2export$eta
Alpha.P <-  items2export$alpha
Deviance.P <- items2export$deviance
Lambda.P <- items2export$lambda
rm (items2export)


load ("objects_2MPE_Mexico.Rdata")
Theta.C <-  items2export$theta
Delta.C <-  items2export$delta
Beta.C  <-  items2export$beta
Eta.C   <-  items2export$eta
Alpha.C <-  items2export$alpha
Deviance.C <- items2export$deviance
Lambda.C <- items2export$lambda
voteDis <- items2export$voteDis
voteAgr <- items2export$voteAgr
voteUnc <- items2export$voteUnc
rm (items2export)

mean (Deviance.cjr); mean (Deviance.P); mean (Deviance.C)


load ("mexico_rollcall_long_format_mean.RData")

# PLOT ideal points
party.color <- as.numeric (by (good.comp.party$party, good.comp.party$final.leg.no, unique))
party.color <- car::recode (party.color, "1='green';
                            2='blue';
                            3='red';
                            4='yellow';
                            5='orange';
                            6='black';
                            7='purple';
                            8='gray'")



# Gather lambda

niceNames <- c("PRI pressure","PAN pressure","PRD pressure","Pressure x PR", "PR"
               , "Pressure x Poverty", "Poverty"
               , "Pressure x Year","Year"
               , "Pressure x Victory Margin", "Victory Margin"
               , "Pressure x State Legislators", "State Legislators"
               , "Pressure x Cabinet Posts", "Cabinet Posts"
               , "Pressure x Presidential Position", "Presidential Position"
               , "Pressure x Committee Chair",  "Committee Chair"
               , "Pressure x PCB Request","PCB Request")

colnames (Lambda.C) <- colnames (Lambda.P) <- niceNames


# Estimate prob that PAN pressure exerts largest effect
set.seed (1971)
PRI.press.coef <- sample (Lambda.C[,grep("PRI pressure",colnames(Lambda.C))], 100, replace=TRUE)
PAN.press.coef <- sample (Lambda.C[,grep("PAN pressure",colnames(Lambda.C))], 100, replace=TRUE)
PRD.press.coef <- sample (Lambda.C[,grep("PRD pressure",colnames(Lambda.C))], 100, replace=TRUE)
PAN.largest <- ifelse (PAN.press.coef > PRI.press.coef & PAN.press.coef > PRD.press.coef, 1, 0)

# Probability that PAN exerts more pressure than other parties
length (PAN.largest[PAN.largest==1])/length (PAN.largest) # prob that PRI effect is larger than PAN effect

# Collect additional information
PRI.theta <- round (median(colMeans(Theta.C)[party.color=="green"]), 1)
PRI.delta <- round (median(colMeans(Delta.C)[party.color=="green"]), 1)

PAN.theta <- round (median(colMeans(Theta.C)[party.color=="blue"]), 1)
PAN.delta <- round (median(colMeans(Delta.C)[party.color=="blue"]), 1)

PRD.theta <- round (median(colMeans(Theta.C)[party.color=="red"]), 1)
PRD.delta <- round (median(colMeans(Delta.C)[party.color=="red"]), 1)

# Print assumed party position
print (PRI.theta); print (PRI.delta)
print (PAN.theta); print (PAN.delta)
print (PRD.theta); print (PRD.delta)

# Legislators ideal points (point predictions)
PRI.theta.legs <- colMeans(Theta.C)[party.color=="green"]
PRI.delta.legs <- colMeans(Delta.C)[party.color=="green"]

PAN.theta.legs <- colMeans(Theta.C)[party.color=="blue"]
PAN.delta.legs <- colMeans(Delta.C)[party.color=="blue"]

PRD.theta.legs <- colMeans(Theta.C)[party.color=="red"]
PRD.delta.legs <- colMeans(Delta.C)[party.color=="red"]

# Fix hypothetical bill position (discrimination and difficulty)
Beta <- 0.5            # we keep item parameters fixed
Eta  <- -0.5
Alpha <- 0

y.star.PRI.party <- Beta*PRI.theta + Eta*PRI.delta - Alpha
y.star.PAN.party <- Beta*PAN.theta + Eta*PAN.delta - Alpha
y.star.PRD.party <- Beta*PRD.theta + Eta*PRD.delta - Alpha

# Print probability that party votes Aye
pnorm (y.star.PRI.party); round (pnorm(y.star.PRI.party)) # votes in favor
pnorm (y.star.PAN.party); round (pnorm(y.star.PAN.party)) # votes in favor
pnorm (y.star.PRD.party); round (pnorm(y.star.PRD.party)) # votes against

# Fix hypothetical directional pressure scores
pressure <- c(0.5,1.5) # two assumed pressure scores (low and high)

# PRI legislators' party-pressure propensity under low positive and high positive pressure?
y.star.PRI <- matrix (NA, nrow = length(PRI.theta.legs), ncol=2)
for (i in 1:length(PRI.theta.legs)){
  for (j in 1:2){
    y.star.PRI[i,j] <- Beta*PRI.theta.legs[i] + Eta*PRI.delta.legs[i] - Alpha + 
      mean(Lambda.P[,grep("PRI pressure",colnames(Lambda.P))])*pressure[j]
  }
}

# PAN legislators' party-pressure propensity under low positive and high positive pressure?
y.star.PAN <- matrix (NA, nrow = length(PAN.theta.legs), ncol=2)
for (i in 1:length(PAN.theta.legs)){
  for (j in 1:2){
    y.star.PAN[i,j] <- Beta*PAN.theta.legs[i] + Eta*PAN.delta.legs[i] - Alpha + 
      mean(Lambda.P[,grep("PAN pressure",colnames(Lambda.P))])*pressure[j] # Negative pressure
  }
}

# PRI legislators' party-pressure propensity under low positive and high positive pressure?
y.star.PRD <- matrix (NA, nrow = length(PRD.theta.legs), ncol=2)
for (i in 1:length(PRD.theta.legs)){
  for (j in 1:2){
    y.star.PRD[i,j] <- Beta*PRD.theta.legs[i] + Eta*PRD.delta.legs[i] - Alpha + 
      mean(Lambda.P[,grep("PRD pressure",colnames(Lambda.P))])*pressure[j]*(-1) # Negative pressure
  }
}


# Now introduce some uncertainty
PRI.prob.aye.voting <- pnorm (y.star.PRI)
PAN.prob.aye.voting <- pnorm (y.star.PAN)
PRD.prob.aye.voting <- pnorm (y.star.PRD)

vote.cond <- function (x) {
  vote <- rbinom (1, size = 1, prob=x)
  return (vote)
}

abs.cond <- function (vote, party=1, agree=0.5, disagree=0.5) {
  if (vote==party) {
    abs <- rbinom (1, size = 1, prob=agree) * (-1) + 1
  } else {
    abs <- rbinom (1, size = 1, prob=disagree) * (-1) + 1
  }
  return (abs)
}

set.seed (1971)

PRI.total.no.yes <- PRI.total.no.abs <- PRI.total.no.no <- c()
for (j in 1:1000){
  vote <- apply (PRI.prob.aye.voting, c(1,2), vote.cond)
  abs  <- apply (vote, c(1,2), abs.cond, agree=mean(voteAgr), disagree=mean(voteDis))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- colSums (vote * obs)
  no.abs <- colSums (abs) 
  no.no  <- colSums (vote.nay * obs)
  PRI.total.no.yes <- rbind (PRI.total.no.yes, no.yes)
  PRI.total.no.abs <- rbind (PRI.total.no.abs, no.abs)
  PRI.total.no.no  <- rbind (PRI.total.no.no, no.no)
}


PAN.total.no.yes <- PAN.total.no.abs <- PAN.total.no.no <-c()
for (j in 1:1000){
  vote <- apply (PAN.prob.aye.voting, c(1,2), vote.cond)
  abs  <- apply (vote, c(1,2), abs.cond, agree=mean(voteAgr), disagree=mean(voteDis))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- colSums (vote * obs)
  no.abs <- colSums (abs) 
  no.no  <- colSums (vote.nay * obs)
  PAN.total.no.yes <- rbind (PAN.total.no.yes, no.yes)
  PAN.total.no.abs <- rbind (PAN.total.no.abs, no.abs)
  PAN.total.no.no  <- rbind (PAN.total.no.no, no.no)
}


PRD.total.no.yes <- PRD.total.no.abs <- PRD.total.no.no <-c()
for (j in 1:1000){
  vote <- apply (PRD.prob.aye.voting, c(1,2), vote.cond)
  abs  <- apply (vote, c(1,2), abs.cond, party=0, agree=mean(voteAgr), disagree=mean(voteDis))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- colSums (vote * obs)
  no.abs <- colSums (abs) 
  no.no  <- colSums (vote.nay * obs)
  PRD.total.no.yes <- rbind (PRD.total.no.yes, no.yes)
  PRD.total.no.abs <- rbind (PRD.total.no.abs, no.abs)
  PRD.total.no.no  <- rbind (PRD.total.no.no, no.no)
}

PRI.aye <- apply (PRI.total.no.yes, 2, quantile, prob=0.5)
PRI.nay <- apply (PRI.total.no.no,  2, quantile, prob=0.5)
PRI.abs <- apply (PRI.total.no.abs,  2, quantile, prob=0.5)

PAN.aye <- apply (PAN.total.no.yes, 2, quantile, prob=0.5)
PAN.nay <- apply (PAN.total.no.no,  2, quantile, prob=0.5)
PAN.abs <- apply (PAN.total.no.abs,  2, quantile, prob=0.5)

PRD.aye <- apply (PRD.total.no.yes, 2, quantile, prob=0.5)
PRD.nay <- apply (PRD.total.no.no,  2, quantile, prob=0.5)
PRD.abs <- apply (PRD.total.no.abs,  2, quantile, prob=0.5)

Table.vote.expectations <- rbind (round (c(PRI.aye[1], PRI.nay[1], PRI.abs[1], PRI.aye[2], PRI.nay[2], PRI.abs[2])),
                                  round (c(PAN.aye[1], PAN.nay[1], PAN.abs[1], PAN.aye[2], PAN.nay[2], PAN.abs[2])),
                                  round (c(PRD.aye[1], PRD.nay[1], PRD.abs[1], PRD.aye[2], PRD.nay[2], PRD.abs[2])))
colnames (Table.vote.expectations) <- c("Lo.aye", "Lo.nay", "Lo.abs", "Hi.aye", "Hi.nay", "Hi.abs")
rownames (Table.vote.expectations) <- c("PRI", "PAN", "PRD")

# Table 6 (main paper)
print (Table.vote.expectations)

########################################################################
########################################################################
### Brazil, probability that parties will vote on a bill
########################################################################
########################################################################

rm (list=ls())

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# Set graph path
graphPath <- getwd()


# Create functions used later
vote.cond <- function (x) {
  vote <- rbinom (1, size = 1, prob=x)
  return (vote)
}

abs.cond <- function (vote, party=1, agree=0.5, disagree=0.5) {
  if (vote==party) {
    abs <- rbinom (1, size = 1, prob=agree) * (-1) + 1
  } else {
    abs <- rbinom (1, size = 1, prob=disagree) * (-1) + 1
  }
  return (abs)
}

load ("objects_2MPE_Brazil.Rdata")

load (file="brazil_rollcall_long_format_mean.Rdata")  

new.party.color <- as.numeric (unlist( by (molten.party.member$party.no, molten.party.member$final.leg.no, unique)))
new.party.color <- car::recode(new.party.color, " 1='red'
                               ; 2='purple'
                               ; 3='orange'
                               ; 4='yellow'
                               ; 5='black'
                               ; 6='blue'
                               ; 7='gray'")


niceNames <- c("PT directional pressure","PSDB directional pressure"
               , "DEM directional pressure","PMDB directional pressure"
               , "PR directional pressure","PP directional pressure"
               , "Directional pressure x Seniority","Seniority"
               , "Directional pressure x District income","District income"
               , "Directional pressure x Commmittee chair", "Committee chair"
               , "Directional pressure x Percentile rank", "Percentile rank"
               , "Directional pressure x District magnitude","District magnitude"
               , "Directional pressure x Year","Year"
               , "Directional pressure x Procedural vote","Procedural vote"
               , "Directional pressure x Cabinet post", "Cabinet post"
               , "Directional pressure x Prez position", "Prez position"
               , "Directional pressure x Party agreements","Party agreements")

Theta.Brz.C <-  items2export$theta
Delta.Brz.C <-  items2export$delta
Beta.Brz.C <- items2export$beta
Eta.Brz.C <- items2export$eta
Alpha.Brz.C <- items2export$alpha
Lambda.Brz.C <- items2export$lambda
colnames (Lambda.Brz.C) <- niceNames
probVoteDisagree  <- items2export$voteDis
probVoteAgree     <- items2export$voteAgr
probVoteUncertain <- items2export$voteUnc
rm (items2export)

# Estimate prob that PT pressure is largest
set.seed (1971)
PT.press.coef   <- sample (Lambda.Brz.C[,1], 100, replace=TRUE)
PSDB.press.coef <- sample (Lambda.Brz.C[,2], 100, replace=TRUE)
DEM.press.coef  <- sample (Lambda.Brz.C[,3], 100, replace=TRUE)
PMDB.press.coef <- sample (Lambda.Brz.C[,4], 100, replace=TRUE)
PR.press.coef   <- sample (Lambda.Brz.C[,5], 100, replace=TRUE)
PP.press.coef   <- sample (Lambda.Brz.C[,6], 100, replace=TRUE)

PT.largest <- ifelse (PT.press.coef > PSDB.press.coef & PT.press.coef > PMDB.press.coef
                      & PT.press.coef > DEM.press.coef & PT.press.coef > PR.press.coef
                      & PT.press.coef > PP.press.coef, 1, 0)

# Probability that PT exerts more pressure than other parties
length (PT.largest[PT.largest==1])/length (PT.largest) # prob that PSDB effect is larger than PT effect

# PT median
PT.theta <- median(colMeans(Theta.Brz.C)[new.party.color=="red"])
PT.delta <- median(colMeans(Delta.Brz.C)[new.party.color=="red"])

# PT eccentric
PT.theta.leg <- median(colMeans(Theta.Brz.C)[new.party.color=="red"]) + 1.95*sd(colMeans(Theta.Brz.C)[new.party.color=="red"])
PT.delta.leg <- median(colMeans(Delta.Brz.C)[new.party.color=="red"]) - 1.95*sd(colMeans(Delta.Brz.C)[new.party.color=="red"])

# Bill characteristics
Beta=0.5; Eta=-0.5; Alpha=-0.5

# propensities to vote
y.star.PT.party <- Beta*PT.theta + Eta*PT.delta - Alpha
pnorm (y.star.PT.party); round (pnorm(y.star.PT.party), 2) # votes in favor

y.star.PT.leg <- Beta*PT.theta.leg + Eta*PT.delta.leg - Alpha
pnorm (y.star.PT.leg); round (pnorm(y.star.PT.leg), 2) # votes in favor

# Add predictors
seniority <- c(0,10)
committee <- c(0,1)
pref.rank <- c(1,0)
dist.magn <- c(2.08,4.25)
year      <- c(1,4)
proc.vote <- c(0,1)
party.agr <- c(0,1)
pressure  <- c(0,-1.5)

# Predictors for pressure on specific legislators
preds <- c("PT directional pressure",
           "Directional pressure x Seniority",
           "Seniority",
           "Directional pressure x Commmittee chair",
           "Committee chair",
           "Directional pressure x Percentile rank",
           "Percentile rank",
           "Directional pressure x District magnitude",
           "District magnitude",
           "Directional pressure x Year",
           "Year")

pred.point.est <- colMeans (Lambda.Brz.C[,preds])

# PT, no seniority, no committee chair, low rank, low magnitude
press <- pressure[1]
leg1.lo <- c(1, press*seniority[1], seniority[1], press*committee[1], committee[1],
             press*pref.rank[1], pref.rank[1], press*dist.magn[1], dist.magn[1],
             press*year[1], year[1])

leg2.lo <- c(1, press*seniority[2], seniority[2], press*committee[2], committee[2],
             press*pref.rank[2], pref.rank[2], press*dist.magn[2], dist.magn[2],
             press*year[1], year[1])

press <- pressure[2]
leg1.hi <- c(1, press*seniority[1], seniority[1], press*committee[1], committee[1],
             press*pref.rank[1], pref.rank[1], press*dist.magn[1], dist.magn[1],
             press*year[1], year[1])

leg2.hi <- c(1, press*seniority[2], seniority[2], press*committee[2], committee[2],
             press*pref.rank[2], pref.rank[2], press*dist.magn[2], dist.magn[2],
             press*year[1], year[1])


prop.leg1.lo <- y.star.PT.leg + pred.point.est %*% leg1.lo
prop.leg2.lo <- y.star.PT.leg + pred.point.est %*% leg2.lo
prop.leg1.hi <- y.star.PT.leg + pred.point.est %*% leg1.hi
prop.leg2.hi <- y.star.PT.leg + pred.point.est %*% leg2.hi


leg1.lo.no.yes <- leg1.lo.no.abs <- leg1.lo.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.leg1.lo))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  leg1.lo.no.yes <- c(leg1.lo.no.yes, no.yes)
  leg1.lo.no.abs <- c(leg1.lo.no.abs, no.abs)
  leg1.lo.no.no  <- c(leg1.lo.no.no, no.no)
}


leg1.hi.no.yes <- leg1.hi.no.abs <- leg1.hi.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.leg1.hi))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  leg1.hi.no.yes <- c(leg1.hi.no.yes, no.yes)
  leg1.hi.no.abs <- c(leg1.hi.no.abs, no.abs)
  leg1.hi.no.no  <- c(leg1.hi.no.no, no.no)
}


leg2.lo.no.yes <- leg2.lo.no.abs <- leg2.lo.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.leg2.lo))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  leg2.lo.no.yes <- c(leg2.lo.no.yes, no.yes)
  leg2.lo.no.abs <- c(leg2.lo.no.abs, no.abs)
  leg2.lo.no.no  <- c(leg2.lo.no.no, no.no)
}


leg2.hi.no.yes <- leg2.hi.no.abs <- leg2.hi.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.leg2.hi))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  leg2.hi.no.yes <- c(leg2.hi.no.yes, no.yes)
  leg2.hi.no.abs <- c(leg2.hi.no.abs, no.abs)
  leg2.hi.no.no  <- c(leg2.hi.no.no, no.no)
}

pred.legs <- data.frame (matrix (data=c(sum(leg1.lo.no.yes)/1000, sum(leg1.lo.no.abs)/1000, sum(leg1.lo.no.no)/1000,
                                        sum(leg1.hi.no.yes)/1000, sum(leg1.hi.no.abs)/1000, sum(leg1.hi.no.no)/1000,
                                        sum(leg2.lo.no.yes)/1000, sum(leg2.lo.no.abs)/1000, sum(leg2.lo.no.no)/1000,
                                        sum(leg2.hi.no.yes)/1000, sum(leg2.hi.no.abs)/1000, sum(leg2.hi.no.no)/1000),
                                 ncol = 3, byrow = TRUE))

colnames (pred.legs) <- c("Aye","Abs","Nay")
rownames (pred.legs) <- c("leg1-lo","leg1-hi","leg2-lo","leg2-hi")

# Table 8a (main paper)
print (round (pred.legs, 2))



#############################################
# Predictors for pressure on specific votes #
#############################################
# Add predictors
year  <- c(1,4)
proc  <- c(0,1)
prez  <- c(1,0)
pressure  <- c(0,-1.5)

preds <- c("PT directional pressure",
           "Directional pressure x Year",
           "Year",
           "Directional pressure x Procedural vote",
           "Procedural vote",
           "Directional pressure x Prez position",
           "Prez position"
)

pred.point.est <- colMeans (Lambda.Brz.C[,preds])

# PT average legislator, year 1, no procedural, no president support
press <- pressure[1]
bill1.lo <- c(1, press*year[1], year[1],
              press*proc[1], proc[1],
              press*prez[1], prez[1])

bill2.lo <- c(1, press*year[2], year[2],
              press*proc[2], proc[2],
              press*prez[2], prez[2])

press <- pressure[2]
bill1.hi <- c(1, press*year[1], year[1],
              press*proc[1], proc[1],
              press*prez[1], prez[1])

bill2.hi <- c(1, press*year[2], year[2],
              press*proc[2], proc[2],
              press*prez[2], prez[2])


prop.bill1.lo <- y.star.PT.leg + pred.point.est %*% bill1.lo
prop.bill2.lo <- y.star.PT.leg + pred.point.est %*% bill2.lo
prop.bill1.hi <- y.star.PT.leg + pred.point.est %*% bill1.hi
prop.bill2.hi <- y.star.PT.leg + pred.point.est %*% bill2.hi


bill1.lo.no.yes <- bill1.lo.no.abs <- bill1.lo.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.bill1.lo))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  bill1.lo.no.yes <- c(bill1.lo.no.yes, no.yes)
  bill1.lo.no.abs <- c(bill1.lo.no.abs, no.abs)
  bill1.lo.no.no  <- c(bill1.lo.no.no, no.no)
}


bill1.hi.no.yes <- bill1.hi.no.abs <- bill1.hi.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.bill1.hi))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  bill1.hi.no.yes <- c(bill1.hi.no.yes, no.yes)
  bill1.hi.no.abs <- c(bill1.hi.no.abs, no.abs)
  bill1.hi.no.no  <- c(bill1.hi.no.no, no.no)
}


bill2.lo.no.yes <- bill2.lo.no.abs <- bill2.lo.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.bill2.lo))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  bill2.lo.no.yes <- c(bill2.lo.no.yes, no.yes)
  bill2.lo.no.abs <- c(bill2.lo.no.abs, no.abs)
  bill2.lo.no.no  <- c(bill2.lo.no.no, no.no)
}


bill2.hi.no.yes <- bill2.hi.no.abs <- bill2.hi.no.no <-c()
for (j in 1:1000){
  vote <- vote.cond (pnorm(prop.bill2.hi))
  abs  <- abs.cond (vote, party=0, agree=mean(probVoteAgree), disagree=mean(probVoteDisagree))
  vote.nay <- vote*(-1) + 1
  obs  <- abs*(-1) + 1
  no.yes <- vote * obs
  no.abs <- abs 
  no.no  <- vote.nay * obs
  bill2.hi.no.yes <- c(bill2.hi.no.yes, no.yes)
  bill2.hi.no.abs <- c(bill2.hi.no.abs, no.abs)
  bill2.hi.no.no  <- c(bill2.hi.no.no, no.no)
}

pred.bills <- data.frame (matrix (data=c(sum(bill1.lo.no.yes)/1000, sum(bill1.lo.no.abs)/1000, sum(bill1.lo.no.no)/1000,
                                         sum(bill1.hi.no.yes)/1000, sum(bill1.hi.no.abs)/1000, sum(bill1.hi.no.no)/1000,
                                         sum(bill2.lo.no.yes)/1000, sum(bill2.lo.no.abs)/1000, sum(bill2.lo.no.no)/1000,
                                         sum(bill2.hi.no.yes)/1000, sum(bill2.hi.no.abs)/1000, sum(bill2.hi.no.no)/1000),
                                  ncol = 3, byrow = TRUE))

colnames (pred.bills) <- c("Aye","Abs","Nay")
rownames (pred.bills) <- c("bill1-lo","bill1-hi","bill2-lo","bill2-hi")

# Table 8b (main paper)
print (round (pred.bills, 2))
